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DISCRIMINANT FUNCTIONS WITH COVARIANCE 


By W. G. Cocuran anv C. I. Buss 


North Carolina State College; Connecticut Agricultural Experiment Station and 
Yale University 


1. Summary. This paper discusses the extension of the discriminant func- 
tion to the case where certain variates (called the covariance variates) are known 
to have the same means in all populations. Although such variates have no 
discriminating power by themselves, they may still be utilized in the discriminant 
function. 

The first step is to adjust the discriminators by means of their ‘within-sample’ 
regressions on the covariance variates. The discriminant function is then 
calculated in the usual way from these adjusted variates. The standard tests of 
significance for the discriminant function (e.g. Hotelling’s T’ test) can be ex- 
tended to this case without difficulty. A measure is suggested of the gain in 
information due to covariance and the computations are illustrated by a numeri- 
calexample. The discussion is confined to the case where only a single function 
of the population means is being investigated. 


2. Introduction. Discriminant function analysis is now fairly well advanced 
for the case where there are only two populations. The data consist of a number 
of measurements, called the discriminators, that have been made on each member 
of a random sample from each population. The technique has various uses. 
Fisher [1] used it in seeking a linear function of the measurements that could be 
employed to classify new observations into one or other of the two populations. 
He pointed out [2] that a test of significance of the difference between the two 
samples, developed from his discriminant, was identical with Hotelling’s generali- 
zation of Student’s ¢ test, discovered some years earlier [3]. Mahalanobis’ con- 
cept of the generalized distance between two populations [4] was also found to 
be closely related to the discriminant function. In any of these applications— 
to classification, testing significance, or estimating distance—we may also be 
interested in considering whether certain of the measurements really contribute 
anything to the purpose at hand, and helpful tests of significance are available 
for this purpose. 

Recently the authors encountered a problem in which it seemed advisable to 
combine discriminant function analysis with the analysis of covariance. This 
case occurs whenever, in addition to the discriminators, there is a measurement 
whose mean is known to be the same in both populations. Suppose, for example, 
that the I.Q.’s of each of a sample of students are measured. The sample is 
then divided at random into two groups, each of which subsequently receives a 
different type of training. Measurements made at the end of the period of train- 
ing would be potential discriminators, but in the case of the initial 1.Q.’s we can 
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152 W. G. COCHRAN AND C. I. BLISS 
clearly assume that there is no difference in the means of the populations cor. 
responding to the two groups. 

The initial 1.Q. measurements are of course of no use in themselves in studying 
differences introduced by the training. Nevertheless, if they are correlated with 
the discriminators, they may serve in some way to ‘improve’ the discriminant: 
e.g. to increase the power of Hotelling’s T” test, or to reduce the number of errors 
in classification. This paper discusses the problem of utilizing such measure- 
ments, which will be called covariance variates. The problem is analogous 
to that which is solved by the analysis of covariance. In covariance, as applied 
for instance in a controlled experiment, variates that are unaffected by the 
experimental treatments can be used to provide more accurate estimates of the 
effects of the treatments or to increase the power of the F test of the differences 
among the treatment means. 

The procedure suggested is as follows. First, the multiple regression is ob- 
tained of each discriminator on all the covariance variates. These regressions 
are calculated from the ‘within-sample’ sums of squares and products: that is, 
from the sums of squares and products of deviations of the individual measure- 
ments from their sample means. Each discriminator is then replaced by its 
deviations from the multiple regression, and a new discriminant function is 
calculated in the usual way from these deviations. The extensions of Hotelling’s 
T’ and Mahalanobis’ distance are both obtained from this discriminant, though a 
further adjustment factor is needed for tests of significance. 

This paper is arranged in three parts. Part I presents a numerical example. 
The decision to place the example first was taken because most of the actual 
applications of the discriminant function in the literature appear to have been 
made by persons relatively unfamiliar with the theory of multivariate analysis. 
It is hoped that with the aid of the example readers in this class may be able to 
utilize covariance variates. For the same reason, the calculations have been 
presented as far as possible in terms of the operations of ordinary multiple re- 
gression, rather than in the form in which they first emerge from the theory. 
Actually, various equivalent methods of calculation are available, and it is not 
claimed that our method is necessarily the best. A mathematical statistician 
may prefer to follow the computing methods which come directly from theory 
(Part II, section 13). 

The example is more complex in structure than the two-sample case. The 
data constitute a two-way classification, in which the row means are nuisance 
parameters, being of no interest, while only a single linear function of the column 
means is of interest. It is well known that the ordinary ¢ test can be applied 
not only to the difference between two sample means, but to any linear function 
of a number of sample means in data that are quite complex. Discriminant 
function technique can be extended in the same way, and readers familiar with 
the analysis of variance should find no great difficulty in making the appropriate 
extension to such data. 

Part II presents the theory. The reader who is primarily interested in theory 
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should read Part II before Part I. Since the approaches used by Mahalanobis, 
Hotelling and Fisher all converge, we have chosen that of Mahalanobis, mainly 
because the extension of his techniques to include covariance variates seems 
straightforward. Maximum likelihood estimation of the generalized distance 
js presented in full for the two-population case. The frequency distribution of 
the estimated distance and the extension of the T° test are worked out. An 
attempt is also made to obtain a quantity that will measure what has been gained 
by the use of covariance. 

In order to illustrate how the theory applies with other types of data, the 
mathematical model is given for the row by column classification that occurs in 
the example. The major results for this model are indicated, though without 
proof. 

In Part III it is shown that the computational methods used in the example 
are equivalent to those developed by theory. While this can easily be verified 
in a particular case, it is not intuitively obvious. 


PART I NumericaLt ExaMPLE 


3. Description. The data form part of an experiment on the assay of insulin 
of which other parts have been published [5]. Twelve rabbits were used. 
Each rabbit received in succession four doses of insulin, equally spaced on a log. 
scale. An interval of eight days or more elapsed between successive doses, and 
the order in which the doses were given to any rabbit was determined by random- 
ization. Thus the experiment is of the ‘randomized blocks’ type, where each 
rabbit constitutes a block and there are 12 blocks with 4 treatments each. 

The effect of insulin is usually measured by some function of the blood sugar 
of the rabbit in periodic bleedings after injection of the insulin. The blood sugar 
was measured for each rabbit at 1, 2, 3, 4, and 5 hours after injection, and also 
before injection. In order to simplify the arithmetic, only the initial blood sugar 
and the blood sugars at 3 and 4 hours after injection will be considered here. 
These data are shown for the first three rabbits (with totals for all 12 rabbit=) 
in Table I. 

Let Xiwz be a typical observation of blood sugar, where 7 = 3, 4 stands for the 
hour after injection, w for the rabbit and d for the dose. The mathematical 
model to be used is as follows. 


(1) Liwe = bi + pPiw + Yis + Bio(Xows — Xo--) + Civs - 


The parameters y; , pi» and yi, represent the true mean and the effects of rabbit 
and log dose respectively. The quantity 2ow: is the initial blood sugar for the 
rabbit w before the test at dose z, while z.. is the average initial blood sugar over 
the whole experiment. The blood sugar at 7 hours has been found experimentally 
to be correlated with the corresponding initial blood sugar, and the relationship 
is represented here as a linear regression, with fi as the regression coefficient. 
The residuals e;,, are assumed to follow a multivariate (in this case bivariate) 
normal distribution, with zero means. The covariance between @ju2 and € jos 
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is taken as o;;.0. The model is the standard one for the ordinary analysis of 
covariance, except that we have two measures of the effect of insulin, z; and a, . 

One additional assumption was made. For all post-injection readings, the 
blood sugar seemed linearly related to the log dose t,. Since this result has been 
found in other experiments, we assumed that 


= bjt, 
where 6; is the regression coefficient of blood sugar on log dose. 
4. Object of the analysis. Our object was to find the linear combination of 


the three blood sugar readings that would measure best the effect of the insulin, 
Because of the linearity of the regression on log dose, the effect of insulin on each 








TABLE 1 
Sample of original data on blood sugar levels in insulin experiment 
| Log ion 
— panes blood omape zo| Three hoursz; | Four hours 24 
32 | aT | 62 | a7 “32 AT | 62 | 77 | 32. a7 | 62) 7 
1 75 94.107, 94 95 76 67| 56 96 95 115 91 
2 91 86 83 93 98 90 77, 69 104 87 90 89 
3 97 99 90 91 84 76 59 48 93 102 85 90 
NE a iicvnsa was 1065 1074 11211070 932 872 731 5911098 1026 970 847 





*12 vabbite. 


x; is known completely if the slope 6; is known. It seems reasonable to choose 
the linear compound of the x;’s which will give the maximum ratio when its 
estimated regression on log dose is divided by the estimated standard error of 
this regression. We now consider how to obtain this maximum. The argument 
given below is not intended to prove the v ny of the method, for which refer- 
ence should be made to Part IT. 

The true regression of the original blood sugar 2) on log dose is known to be 
zero. Hence, it is clear that the variate xo is useful only in so far as it enables 
us to obtain more accurate estimates of 6; and 6,. For this purpose we need to 
estimate the effect of xo upon 23 and 2, , the blood sugar readings at 3 and 4 hours, 
independently of dose of insulin or of differences between rabbits. From the 
standard theory of covariance the best estimate is the regression coefficient 
bio = Eio/Eoo , where E denotes a sum of squares or products calculated from 
the error line in the analysis of covariance; that is from the sums of squares and 
products of deviations of the x; from the fitted regression on row and column 
parameters. 
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The regression of the blood sugar at each hour on the log dose of insulin is 
calculated from totals adjusted for the regression on 2). Since the 4 successive 
log doses (z = 1, 2, 3, 4) are spaced equally, they may be replaced in the com- 
putation by the coded doses —3, —1, +1, and +3. If we let 7;, be the total 
blood sugar, summed over 12 rabbits, at the 7th hour with dose z, the following 
result is well known for the analysis of covariance. The best estimate of 
6(2 = 3, 4) is 

((-37 a — Tie t+ Tis + 387i) — bio (—3Tn — Toe + To: + 37'04))/240. 
The divisor, 240, is 12(37 + 17+ 1° + 3°). The expression may be written 
di _ (di — bods) 
240 240 : 
where 
dj = —3Ta — Te t+ Ts + 3TH. 


A linear combination is formed from d3; and d;, the numerators in the best 
estimates of 5; and 6, , by means of the coefficients Z; and L,. lL; and L, are 
computed so as to maximize the ratio of 


d,; = L3d3 + Lid 


to its estimated standard error. 
From the definition of d; , this requires a discriminant of the form 


I = L3(X3we a b30Xowe) + Ly ( ave — b.oLows), 


where each 2,2 is measured from its mean. 
We require next the estimated standard error of d;. This depends, in turn 
. , / . ° . ° 
upon the variances of d; and d, and their covariance. As usual in the analysis, 
of variance we have 


2 
(5) V (ds) = V(d;) + dj V(b30) = 233.0 (240 + é) ° 
00 
The residual variance o3;.9 is estimated from the sums of squares and products 
in the error row of the analysis of covariance as 
83.0 = Ey3.0/n = (Ess — E%o/Ev)/n, 


where 7 is the degrees of freedom in each E;; diminished by one. Similar methods 
lead to the variance of d; and to the covariance of d; and d;. It follows that the 
true variance of d; may be written 


(6) V (d:) a L; 033.0 + 2L3L, o31.0 + Li 044.0 » 


2 
=) in equation (5) is omitted since it does not involve 


where the factor (210 ae 
Eo 
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the L’s. Similarly, the estimated variance of d; , apart from constant factors, 
may be written as 


(7) L3E339 + 2LsL:Ey9 + LiEuo. 
The quantity to be maximized is therefore 
nds + Iyds) 
V1 Bisg + QL, Bus + LBs 





Formally, this is the same type of quantity that is maximized in ordinary analysis 
with the discriminant function. Differentiation with respect to the L’s leads 
to the equations (after omission of another constant factor) 


(8) Ex; oL3 + Es; ols — d; ’ E3; 9L3 + Ess ols = dj . 


The objective of the computation, therefore, is to obtain discriminant coefficients 
having the same ratio to each other as L; and Ly, in equations (8). As will be 
shown in the next section, this can be accomplished in practice more conveniently 
by substituting an alternative set of three simultaneous equations for the two 
in equations (8). 


5. Calculations. The first step is to form the sums of squares and products 
in the analysis of covariance. With 12 rabbits and 4 doses, the conventional 
breakdown of each total sum is into components for rabbits (11 d.f.), doses 
(3 d.f.) and rabbits < doses (33 d.f.). Because of the assumed linear regression 
on log dose, the sum of squares for doses was further divided into two com- 
ponents. The first (1 d.f.) is the contribution due to this regression. For z;, 
the sum of squares due to regression is d;/240, or in the case of x; , (1164)°/240, 
or 5645. The remaining component, (2 d.f.) is called the curvature, since it 
measures the effect of deviations from the linear regression. ‘The sum of squares 
for curvature is found by subtraction. 

The following points may be noted. (i) For both 2; and 2, , the F ratio of the 
curvature mean square to the rabbits * doses mean square will be found to be 
less than 1, so that the data do not suggest rejection of the hypothesis of a linear 
regression on log dose. (ii) The F ratios of the regression mean squares to the 
rabbits < doses mean squares are highly significant, being 57.8 for x; and 28.7 
for z,;. This indicates, incidentally, that the three-hour reading may be a more 
responsive measure of the effect of insulin than the four-hour reading. (iii) With 
2), the F ratio does not approach significance for either the regression or the 
curvature, as is to be expected. 

A consequence of the assumption of linear regression on log dose is that the 
curvature mean squares and products are estimates of the same quantities as 
the rabbits < doses mean squares and products. Consequently, the lines for 
curvature and rabbits X doses in Table 2 could be added to give 35 d.f. for the 
‘error’ sums of squares or products, F3; , etc. We decided, however, to estimate 
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the error only from the 33 d.f. for rabbits X doses. This was done because it 
seemed to facilitate a test of the curvature of the final discriminant J. (This 
test will not be reported here.) 


The L’s could now be obtained from equations (8). In this case the first 
equation would contain the terms 


Es3.9 = 3223 — (1259)?/2351; E349 = 1200 — (1259) (1340) /2351; 
Ei = th — bod = ~1804 — (Fe) 62, 
leading to the simultaneous equations 
2548.8 Ll; + 482.4 Ly, = —1197.2 
482.4 L; + 2373.2 L, — 844.3, 
which give L;/L; = —.41848/—.27070 = 1.5459. 


TABLE 2 
Sums of squares and products 





Component <¢ x? xz | «Xots ots 3% 


Between rabbits.............| ~~ 9376 11165 | 1952 | | 

Between doses | 58 >) 2810 —247 | | 3981 
(Reg. on log dose........./ | > 5645 2727 —301 3924 
Curvature 83 54 | 

Rabbits X doses | | 3137 | 1259 | 


17112 | 2964 | 3719 - 


Instead of using these equations, we propose to solve alternatively the set of 
three equations 


Sool + Sosks + Sods = dy 
(9) SaoLo + S3sli3 + Ssil, = ds 
S:oLo + Sal; + Sal, — d;, 


where each S;; (¢ = 0, 3, 4) is the sum of squares or products formed by adding 
the error line in the analysis of variance to the line for regression on log dose. 
Thus S;; has 34d.f. The ratio of L; to L,; , as found from equations (9), is exactly 
the same as that found from the original equations (8), as is proved in section 18. 
Further, the solution of the new equations seems to be more useful for performing 
tests of significance, as will appear in following sections. 

Accordingly, the first step after forming the analysis of variance is to set up 
the three equations (9). 
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The equations are solved by means of the inverse matrix. The values of d; 
on the right side of the equations are replaced successively by 1, 0, 0 by 0,.1, 9 
and by 0, 0, 1 to obtain the three sets of values for Lo , L3 and L,. These results 
are given in the first three columns of Table 4 and are designated as ¢;; . 

The L’s follow from the c;; by the usual rule for regressions. For example, 


Ls; = {(.003209)(62) + (.227781)(—1164) + (—.199655)(—809)}-10" = 
— .103417 


TABLE 3 
Equations for determining L3; and Ly. 


958L>. + 8868L; + 5124L, = —1164 


1131L0 + 5124L; + 58641, — 809 


The composite response or discriminant, adjusted for the covariance variate, is 


now taken as 
I = L3 (« —_ Ev 2) oe Ly («x = E' r») 
Koo Koo 
or 


—.103417 («: —— ») — 066883 @ . >) 


9351” 351° 
= .0935032% — .10341773; — .066883z,. 


Note that the value of Ly is not used at this stage and that L; /L, = 1.546 agrees 
with the value found from equations (8). 


TABLE 4 
Inverse matrix (X 10°) and L’s 
(10° ¢;;) di he 
.465408 .003209 — .092568 62 . 100008 
.003209 227781 — .199655 —1164 — .103417 


— .092568 


— .199655 . 362846 —809 — .066883 


A similar method may be followed when there are more discriminators or more 
° e _ ° ° : ° 
covariance variates. With two covariance variates, x) and 2 , for instance, the 
adjusted discriminant would be 


L;3(x3 — by% — b3or0) + Ls(X4 — by% — bio) 


where bz , b39 are the partial regression coefficients of x; on 20 , xo respectively, 
determined from the error line, and similarly for z,. Further, since any linear 
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function of the column (dose) means may be represented as a regression on some 
variate t, , this method may be applied to any linear function of the column means 
in which we are interested, provided that the mathematical model is appropriate. 


6. Test of the regression of the adjusted discriminant on log dose. The 
numerator of the regression of J on the coded doses is 


d; = I3(d; ad bsodo) + Li(dy = baodo). 


Since the regressions of x; and 2, on the coded doses were both significant, it 
may be confidently expected that the regression of J will also be significant. The 
test of significance will, however, be given in case it may be useful in other appli- 
cations. For those who are familiar with multiple regression, the test is perhaps 
most easily made by means of a device due to Fisher [2]. 

Construct a dummy variate y,,. such that y,. is always equal to ¢, , or in our 
case to the coded doses. That is, y takes the value —3 for all observations at 

















TABLE 5 
Analysis of y* and yx, 

d. f ¥ yXi 
eed islitaiaio Sie aiden laid ae 
SD 8a oad nates ite al biuieedinn es ark he 4 11 0 0 
i ol se Lieriedthehasuded 3 240 d, 

Regression on log dose. ................. ; 1 240 d; 
ID. vv sadac cess sated eines Sietiale 2 0 0 
Rabbits X doses = error... eye 33 0 0 
Sum = Error plus reg. on log dose......... 34 | 240 d; 
tet ciedi ot b> li ot bie eae aren Blea dee ee 47 240 d; 


the lowest dose level, and —1, +1, and +3 respectively for all observations at 
the successive higher dosage levels. We shall show that equations (9) solved 
in finding the L’s are formally the same as a set of normal equations for the linear 
regression of y on 2%, 73, and %. 

The following analysis for y* and yx; may easily be verified. 

It will be noted that the sum of products of y and 2; in the sum line is d,. 
Further, S;; is the sum of products of x; and x; for this line. It follows that the 
normal equations for the regression of y on the x’s, as calculated from the “‘sum”’ 
line, are 


SioLo + SisL3 + Sisk = dj (« = 0, 3, 4). 


These are just the equations solved in obtaining the L’s. Consequently, Ls 
and L, are the partial regression coefficients of y on z; and 2,. A test of the null 
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hypothesis that the true values of L; and Ly are both zero can be made by the 
standard method for multiple regression, as will be shown later from theory, 
This test is equivalent to a test of the hypothesis that the true value of d, is zero, 

To apply the test, we require three items in the analysis of variance of y. 
First, the total sum of squares for the Sum line, already seen to be 240 (Table 5), 
Second, the reduction due to a regression on all variates (covariance variates 
plus discriminators). By the usual rules for regression, this is (from Table 4) 


+ (—.066883)(— 809) = 180.69, 


Finally, we need the reduction due to a regression on the variates that are not 
being tested, i.e. on the covariance variates alone. From Table 4, the reduction 
































TABLE 6 
Analysis of variance of dummy variate y 














se line inkl oe eo 
Reduction to regression on covariance variates.| 1 1.62 
Additional reduction to regression on dis- 

a Sol aa eit ela lesd toe 2 179.07 | 89.54 


eee 31 59.31 | 1.913 


Total (from Sum line).................... 34 240.00 


in this case is simply d5/ Soo or (62)"/2367, or 1.62. The difference, 180.69 — 1.62, 
represents the reduction due to the regression of y on L; and Jy, after fitting x. 
The resulting analysis is given below, the degrees of freedom being apportioned 
by the usual rules. 

The F ratio, 89.54/1.913, or 46.80, with 2 and 31 df., is used to test the null 
hypothesis that the adjusted discriminant has no real regression on log dose. 





7. Test of particular discriminators. Another useful test is that of the null 
hypothesis that a particular discriminator, or group of discriminators, contribute 
nothing to the adjusted discriminant. In other words, this is a test of the null 
hypothesis that the true values of a subset of the L’s are all zero. The test is of 
interest in the present investigation, since it would be useful to know whether 
all five hourly readings of the blood sugar are really helpful. As might be 
expected by analogy with the previous section, the test is made by calculating the 
additional reduction due to the regression of y on the particular subset of the 
L’s in question. . 

The test will be illustrated with respect to L;. One method of making the 
test is to re-solve the normal equations with L, omitted. From this solution 
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the reduction due to a regression of y on 2 and 2; alone is obtained. The addi- 
tional reduction due to a regression on 2, is found by subtraction from 180.69. 

However, the additional reduction can be found directly from the well-known 
regression theorem that it is equal to Li/cu. The c’s have already been found 
in Table 4. The result is (.066883)"/(.000362846), or 12.33. This value is 
tested against the residual error mean square of 1.913, F having 1 and 31 df. 
The contribution is found to be significant. 

In fact, by this process a kind of estimated standard error can be attached to 
each of the L’s for the discriminators, using the formula s\/c;; , where s is the 
residual root mean square. Thus for L; , (—.103417), the ‘standard error’ is 
4/ (1.913) (.000227781), or .0209. It should be stressed that at this point the 
analogy with regression is rather thin. The L’s are not normally distributed, 
nor do the estimated standard errors follow their usual distribution. It is, 
however, correct that if the true value of ZL, is zero, L4/s~/ cx, follows the ¢ distri- 
bution with 31 df. Thus, if omission of some discriminators seems warranted, 


TABLE 7 
Analysis of variance for regression of y on the discriminators 


a. f. S. S. M. S. 
Regression 159.20 | 79.60 


a ia i aka cal et dt al le ta ae onde ad 80.80 | 2.525 


these ¢ ratios are relevant in deciding which variate to eliminate first. Strictly 
speaking, the c’s should be re-calculated after each elimination before deciding 
which other discriminators might also be discarded. 


8. Estimation of the gain due to covariance. ‘The tests given above enable us 
to state whether the discriminators contribute significantly, in the statistical 
sense. It is also of interest to investigate what has been gained by the use of 
the covariance variates. From the practical point of view, the question: ““What 
is the gain from covariance?’ might be re-phrased as: ‘‘If x is ignored, how 
many rabbits must be tested in order to estimate the regression on log dose as 
accurately as it was estimated with the adjusted discriminant for 12 rabbits?” 

The theoretical aspects of the question are discussed in section 16; the calcula- 
tions are described here. The only new quantity needed is the F ratio for the 
regression of y on the discriminators alone. This can be obtained by a new 
solution of the normal equations, this time with the covariance variates omitted. 
With just one covariance variate, it is quicker to use the fact that the additional 
reduction to the regression of y on 2, after fitting z; and x, is Li/co, 
or (.100008)*/(.000465408) or 21.49. Consequently, the reduction due to a 
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regression of y on x3 and 2, alone is 180.69 — 21.49, or 159.20. The F ratio, 
79.60/2.525, is 31.52, whereas the F ratio with covariance is 46.80 (from 


Table 6). The quantity suggested from theory for comparing the two techniques 
is 












(m2—2)F _ 
Ne 


1 





where mz is the number of d.f. in the denominator of F. These values are 
{(30 X 31.52/32) — 1} or 28.55 with no covariance and {(29 X 46.80/31) — 1} 
or 42.78 with covariance. The relative information is estimated as 42.78/28.55, 
or 1.50, so that the use of covariance gives 50 per cent more information. In 
other words about 18 rabbits would be needed if the initial blood sugars were 
ignored. To a slight extent this estimate favors the covariance analysis, since 
it ignores the increased accuracy that would accrue from the extra error df. if 
18 rabbits were used without covariance. 


PART II THrory 


9. Notation. The theory will be given first for the two-population case. 
We suppose that a random sample of size N has been drawn from each popula- 
tion. A typical discriminator is written 2;,,~2 and a typical covariance variate 
Xiva » Where 









t,j = 1,2,- - - p denote discriminators, 
£,y = 1,2,- - - k denote covariance variates, 
w = 1, 2 denotes the population, and 


a = 1,2,--- N denotes the order within the sample. 






The population mean of tive is win , and the corresponding sample mean is jy. 
The difference (ui2 — pir) is denoted by 6; and the corresponding estimated 
difference (xj2. — Xi.) by d;. 


















10. Discriminant functions and generalized distance. Since we propose to 
approach the theory by means of the generalized distance, it may be well to 
review briefly the relation between the discriminant and the generalized distance. 
In the ordinary theory (with no covariance variates) it is assumed that the var- 
jates Zina follow a multivariate normal distribution, and that the covariance 
matrix o;; between Ziwa and Xj.a is the same in both populations. The gener- 
alized distance of Mahalanobis is defined by 


P ** :< , — 
(10) pA’ = y % o” 6; 4; . where (oc?) = (o:3) 7 
i,j=1 
In order to estimate this quantity from the sample, we first calculate the mean 
within-sample covariance s;;, where 





(11) sj = 7 > (Ziwa ps Saw See — Ljw-)/2(N aid 1), 


wel a=! 
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The estimated distance is then taken as 


(12) pD* = i s” did; . 
Apart from a factor N/(N — 1), this is the maximum likelihood estimate. 

In the discriminant function used by Fisher (1), the object is to find a linear 
function Iwe = UM ixive , Where the M; are chosen to maximize the ratio of the 
sum of squares between samples to that within samples in the analysis of variance 
of J. This is equivalent to maximizing the ratio of the difference between the 
two sample means of J to the estimated standard error of this difference. As 
Fisher showed (2), the M; (apart from an arbitrary multiplier) are given by 


p 


M; = 7 s d; 
j=t 
Consequently, the difference between the two sample means of /, the discriminant 
function, is 


P 


P. ° 
> Mid; = D> s*did;. 
dank Sil 
This is exactly the same as pD” in equation (12). Thus the discriminant func- 
tion leads to the estimated distance, and vice versa. 


11. Extension to the present problem. In our case there are (p + /) variates 
(p discriminators, k covariance variates) from which to estimate the distance. 
All variates, 2ive and Xwa , are assumed to follow a multivariate normal distribu- 
tion. The covariance matrix, assumed the same in both populations, now has 
(p + k) rows and columns, and may be denoted by 


Oi CF 
(13) d -( ‘ 
Ory Tk 


For each of the covariance variates, it is known that the population means 
ti, Mee are equal, so that the difference 6; is zero. This is the fact that dis- 
tinguishes the problem from ordinary discriminant function analysis. 

Hence, the generalized distance, as defined from all (p + k) variates contains 
no contribution from terms in 6; and is given by 


(14) (p + kA? = Do oth in 5:4). 
1 


i.j= 


The matrix o(/,+«) is that formed by the first p rows and columns of the inverse 
of A. Note that in general this will not be the same as the matrix o” , which is 
the inverse of o;;. 

In the next section we consider the estimation of this quantity from the sample 
data. By analogy with the previous section, it might be guessed that the 
estimate would be of the form =s¢/,4,) did;. The maximum likelihood estimate 
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turns out to be of this form, except that instead of d; we have d; , the difference 
between the two sample means of the deviations of x; from its ‘within-sample’ 
linear regression on the 2; . 


12. Estimation of the distance. It is known that the generalized distance 
is invariant under non-singular linear transformations of the variates. For 
. . / 
convenience, we replace the 2,2 by variates Xiwa , Where 






















k 

/ 

Liwa = Liwa — 7 Bee (Le wee a Miw)- 
fal 


Thus 2jwa is the deviation of x;,,. from its population linear regression on the 


Ztwa- The population mean of 2. is clearly yi, , and the difference between 
the two population means is therefore 6; . 


. . g . 
The covariance matrix of the twa , L:wa May be written 


(15) 





where o;;.: denotes the covariance matrix of the deviations of the 2; from their 
regressions on the 2,2. It follows that in terms of the transformed variates 
the generalized distance is given by 


p 
(16) ptkhar?= > o*5,5, 
btu 
where o”’* is the inverse of the p X p matrix o4;.; . 
The joint distribution of the 2N observations on each of the tina and tue 
is as follows: 


—} k js tN g +) , 
(27) — | e* ae a Td iwa ditewo 


( 2 ON 


Pp os , _ 
exp ’ Z. : . ¥ oo (tia aes Miw) (2 jwa — Ujw) + 


wel a=l t,j=1 


bol 


2 ON E \ 
De De ota — Mew) (Sma — tw) | ; 
ul Gul fagut J 
where o*” is the inverse of the k X k matrix og, . 

We now proceed to estimate A’ in equation (16) by maximum likelihood. 
For this, we obviously need the sample estimates of the o”* and the 6; , and it 
will appear presently that the sample estimates of the 6;; are also required. 
However, it happens that the o*” and the y;, are not needed. Hence the rele- 
vant part of the likelihood function is 





D 


‘ ij: l << . Sfety 7 sk 
(17) L=N log | ot | ~ yp ie Zz a” *( awe o tir) (2 j00 aa jw) 


“ w=1 a=1 i,j=1 
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k 
, 
Viwa = Liwa — ZZ Bit (Ltwe = Mew) ° 
g=1 


Differentiating first with respect to yp; , we obtain 


N Pp m 
(18) > 2 6 (aice — fee) = 0. 


acl j=1 


Except in the case (with probability zero) where our estimate of o”’* turns out 
to be singular, these equations have no solution except 


N 
(19) d (xiwa — Aju) = 0 
for every j, w. Consequently 


a , 
Kjno = L jw 


so that 
A / / ‘ 
6; = bj — An = in. - 2. = d; ro 2. Bit d,. 


This shows that the 8; must also be estimated. Now 


2 2 


OL = - oL Ox ™ x B ij.t ’ 
ae SR ——— = Liwa —~ Mtw joa —~ Miw/)- 

OB iz 2d. 2 OX wa OBiz xz 2 2, " ( : Me Ma; Mi ) 
Once again, unless the estimate of o* is singular, the only solutions of the equa- 
tions formed by equating this quantity to zero are 


2 N 
(20) a , (tive = Mew) (2 jee — Biju) = 0 
for every &, 7. 
Since fj. = Xjo-, the term in yz, vanishes. Substituting for x’ in terms of z 
from (17), we obtain 


2 y 


N » 
: Zz Liwa {ee _ Lijw-) — z Dig Tawa ~ re) = 0 
w=1 a=l g=1 
where b;, stands for the maximum likelihood estimate of 8;,. These equations 
may be written 

k 
(21) x bin Ley = Lie 

n= 


where E denotes a sum of squares or products of deviations from the sample 
means, containing 2(N — 1) degrees of freedom. The equations are therefore 
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the ordinary normal equations for the ‘within-sample’ multiple regression of 
Liwa ON the Liwa. - 
Finally, differentiation of L with respect to the o”* leads to 


2 N 


(22) 2NG4;.: — z. 2. (tive _ Liv) (2 juve — jw:)- 


w=1 a=1 


























This is just the ‘within-samples’ sum of squares or products of the variates 2’. 
On substituting for the x’ in terms of the x and using equations (21), we obtain 


k 
2NGi;.: = E;; = 2. biz Ej: = Ei3-: (say). 
To summarize, the estimated distance is given by means of the equation 


P os AA P ee , 
(p+k)D? = >) 6,8, = 2N Do EY? did;, 


itu ‘jaa 


where E”’* is the inverse of E;;.; and 
k 
dj; = d; = ™™ d:. 
E=1 


This estimate was obtained by assuming all variates jointly normally dis- 
tributed. From the form of the likelihood function (17) it can be seen that the 
M.L. estimate of the distance remains the same under the less restrictive assump- 
tions that the x: are fixed, while the deviations of the ;,,. from their regressions 
on the 2%:,a are jointly normal. 
















13. Computational procedure. An orderly procedure for calculating the 

generalized distance will now be given. From this, the method for computing 

the corresponding discriminant function will be shown. The computations also 

lead to the generalization of Hotelling’s T°. The steps are as follows. 

(i). First form the ‘within-sample’ sums of squares and products of all variates, 

with 2(N — 1) degrees of freedom. ‘These are the quantities denoted by 
E;; ; Ex: ; Eey 

(ii). Invert the matrix E;, , giving E*”. 

(iii). The regression coefficients b;; , estimates of the 6;; , are now obtainable by 

means of the relations 


k 
1 wet 
biz a > E,, E '. 
n=1 
as is clear from the usual matrix solution of equations (21). 


(iv). The sums of squares and products of the deviations of the z; from their 
‘within-sample’ regressions on the x; are now computed from equations (22) 


k 
2NGi;-¢ = Ei;.z _ Ei; aia 2d bi Ex. 







DISCRIMINANT FUNCTIONS 167 


(v). The final step is to invert the matrix E;;.; , giving E**, and to form the 
product 


P - k 
(p+ k)D? =2N >> E** did; where di; = d; — Dd. diz dy. 
gal 


t,j=1 


When there were no covariance variates, the discriminant function J had the 
property that the difference between the two sample means of J was equal to the 
estimated distance (Section 10). This relationship can be preserved when co- 
variance variates are present by defining J so that 


k 
Few = > M; (ie staal pm biz Zea) 
ink funk 
and calculating the weights M; from the equations, 
“ / 
Zz Ej;.. M; = d;. 
For in that case, 


P. -s 
M; on Z E'*§ d;. 


j=1 


Consequently the difference between the two sample means of 7 is 


> Md, = > E** di dj, 


i=l i,j=1 


which (apart from the constant 2N’) is equal to (p + k)D’. 


14. Distribution of the estimated distance. In the ordinary case, with no 
covariance variates, the frequency distribution of the estimated distance has 
been given by several authors, e.g. Hsu [6]. It will be found that in our prob- 
lem the distribution is essentially the same, except that the quantity D’ must be 
multiplied by a new factor and that one set of degrees of freedom entering into 
the result must be changed from (n — p+ 1) to(n —-p—k+1). 

Thus far we have assumed that all variates jointly follow a multivariate nor- 
mal distribution. It is convenient at this stage to regard the covariance variates 
Tiwe aS fixed from sample to sample, and to use the conditional distribution of 
the 2.2, subject to this restriction. It is well known (e.g. Cramér [7, section 
24.6]) that this conditional distribution is the multivariate normal 


on)" [ot TT div 


(23) ee. 
; exp {—} | 2. Z z. 0°  (tiwa —~ Kiw ~~ Viwa) (Ljwa —~ Bjw ~ vi) |b 


wel a=1 t,j=1 


where 


k 
Yiwa = 2 Bit (Zewa Perr Mew). 
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Since the estimated distance is a function of the quantities E;;.; , di, we now 
find the joint distribution of these variates. The joint distribution of the 
sums of squares and products £;;.; is obtained by quoting a slight extension of a 
result due to Bartlett [8], which may be stated as follows. 

Let the variates Xiwa follow the distribution (23) and let 


2 N 


(t) Ei; = Z. ‘ tase 4 Sie Mien i Tin) 


w=l a= 


be a typical ‘within-samples’ sum of squares or products, 
(it) by = BB” 
= 
be the ‘within-samples’ partial regression coefficient of x; on x; , and 
(222) Ei. = Ei; — > bi Ei: 


be the sum of squares or products of deviations from these regressions. Then 
(a) the quantities E;;.: follow the Wishart distribution 


Pp os 
c| Eiz¢ |"? exp \-! ae og} I] dz;;.; 


t,jm1 
with (n — k) d.f., wheren = 2(N — 1), 


(b) thts distribution is independent of that of the bi: , and 
(c) both distributions are independent of that of the means xin. and consequently 
of that of the difference d; = (xi2. — Xi.). 

The result was proved by Bartlett for a sample from a single population. The 
extension to the case of two populations is straightforward and will not be given 
in detail. 

From (b) and (c) it follows that the distribution of the E;;.¢ is independent 
of that of the quantities 


k 
di = d; — Do dig dy. 
f=1 
Further, with the x; variates fixed, the d‘ are linear functions of the zi. with con- 
stant coefficients and hence. follow a multivariate normal distribution, Wilks 


[9]. We now find the means and the covariance matrix of this joint distribution. 
From the joint distribution (23) of the z;,,. , it is easily seen that 


k 
(24) E(d;) = 6 + 2 Big de. 


Also, since by standard regression theory the b,: are unbiased estimates of the 


Biz ’ 


k k 
BY: bie a = > Bie de. 
t=1 t=1 
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Hence, by subtraction, 
(25) E(@di) = &. 


Now 


k k 
Cov (d; d}) = Cov (d; — ae biz d:)(d; — Zz bin d,). 
g=1 n=1 


By (c) the distributions of the d; , b;, are independent, so that there will be no 
contribution from products of the form d;b;,. Hence 


k 
(26) Cov (di d;) = Cov (d; d)) + 2, dz dy Cov (big bj). 
Oh Dd 


Since d; is the difference between the means of two samples of size N, Cov 
(d.d;) is 2 o;;.:/N. The covariance of b;: and bj is more troublesome. Writing 
the expressions for these regression coefficients in terms of the original data, we 
have 


k 
Cov (bigbin) = 2Y EE” Cov (Ex Ej) = 
A,y=l 


k N 


. ' 
rE 
Do EXE" DL rea — Trw-) (ret — Lr.) COV (tive Liat). 
Awl w,zel af 
Since successive observations are assumed independent, the covariance term van- 
ishes unless w = z and a = §, in which case it equals o;;... Thus 


k 
Cov (bigbj,) = oi3-2 2, EME" Ey, = o4;..E™. 
Finally, from (26) 
— 2 - 
(27) Cov (d; d;) = Cij-¢ (< + 2, E" d: dy) = VOiz-k (say). 


Having obtained the distributions of the E;;.; , d;, we may apply Hsu’s result 
[6] for the general distribution of Hotelling’s T’. In our notation, this may be 
stated as follows. 

If the variates di/~/v follow the multivariate normal distribution with means 
8:/+/v and covariance matrix oi;.,, and if the variates E;;.; follow the Wishart 
distribution with (n — k) d.f. and covariance matrix o;;.¢ , the two distributions 
being independent, then 

P 


y= Dd) E* di dj/s, 


#,j=1 


follows the distribution 


© h 1 


28 TE alec tieriinacaatiticia ania aiaaaiaeeiaii 
me pmo h! Blap +h, 3(n —k —pt1)} 


‘ yer + gy tenes dy, 
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P. -* 
T= 4 be oa” §5;8;/v, 


4, jal 


2 k 
oe et fae, n = 2(N ~ 1) 


This distribution is, of course, the distribution of the ratio of two independent 


values of x’, with p and (n — k — p + 1) df. respectively, in the case where the 
numerator is non-central. 


15. Tests of significance. This result leads to the extension of Hotelling’s 
T’ test. For if 6; = 0, (¢ = 1,2,--- p), then 7 is zero and 


Set deal 

t,j=1 
is distributed as vpF/(n — k — p+ 1), with pand (n —k — p+ 1)df. The 
distribution (28) above gives the power function of this test. 

We may also wish to apply a test of this type to a subgroup 2; of the dis- 
criminators (¢ = 1,2,---q <p). Speaking popularly, this is a test of the null 
hypothesis that the above variates x; contribute nothing to the discrimination 
between the two populations, given that the remaining discriminators and the 
covariance variates have already been included.’ To see what is meant more 
precisely, consider the following transformation: 


6 = 2 — D Bar, — 2 Bete, t= 1,2,---4q; 
Zr — 2 Bez:, l q oe l,--*p; 
Xe = 2X, £=1,2,---k, 


where the 6’s are population regression coefficients. Then it is not difficult to 
see that the distance is now given by 


4 o2 , P m- 
(pth = >> os5,+ Do’ 526,, 
t,jm1 lym=q+1 
where co” " is the inverse of the covariance matrix of the deviations of the 2; 
from their regressions on the x; plus the a; , and 
5; = 6; — 2 Bud. 


Consequently if 6; = 0, (i = 1, 2, --- q) the distance is exactly the same as it 
would be if the variates x; were omitted. The test in question is therefore a test 
of the null hypothesis that 6; = 0, (¢ = 1, 2,--- q). 

If both the remaining discriminators x, and the covariance variates 2; are re- 
garded as fixed, the method of proof in the previous section provides an F test 


1 The test is illustrated in section 7. 
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for this hypothesis also. It is found that the sums of squares or products E”’* 
follow a Wishart distribution with (n — k — p + q) df., while the quantities 


Pp k 

dj=da- ym bir dy — D> Dy de 
l=g+l1 §=1 

are normally distributed, with zero means when the null hypothesis is true. This 

leads to the result that 


g ee 
>, Et" did; 


ij=l 
is distributed as v’gF /(n — k — p+ 1), with q and (n — k — p+ 1) df., and 


y = 2 +2zE* d, dz, 

N 
the sum extending over both the covariance variates and the discriminators that 
are not being tested. 


16. Discussion of the gain due to covariance. In this section we attempt to 
construct a measure of the amount that has been gained by the use of the co- 
variance variates. Only a preliminary discussion will be given: a complete dis- 
cussion would be rather lengthy, owing to the many different uses to whichthe 
discriminant function can be put. Perhaps the problem can most easily be seen 
by considering the effect on Hotelling’s generalized T” test of significance. 

The power function of this test, as obtained from equation (28) section 14, de- 
pends on four factors; the level of significance that is chosen, the degrees of free- 
dom »; and nz in the numerator and denominator of F, and the parameter r. 
If the covariance variates were ignored, the usual 7” test could be applied to the 
discriminators alone. In this case we would have 


, . ie ; 
ni = p, ng=n—-ptl, r’ = 420"5,6;/v’, where v’ = 2/N. 
With the covariance variates, we have 


mi = DP, m=n-p-—k+1, T= 430" *5,5;/v, 


b= s + SE" dz dp. 


The first point to note is that 
> 0 5,5; > V 056; 


This is an instance of the general result that the addition of new variates cannot 
. 2 oan . ‘ ° ° 
decrease the value of pA’. To see this, replace the covariance variates by their 
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deviations from their regressions on the discriminators. This transformation 
gives 


Pp *s P ss k * , 
(29) z o? *5, 6; = 7 o 7545; + 2 o*” 5; 5, ’ 
= 


i,jul i,j=l 


where 


Pp 
5: = Og = 2 Beds. 


Since the term on the right of equation (29) is a positive definite quadratic form, 
the result follows. 

Consequently, the first effect of the covariance variates is to make the numer- 
ator of 7 greater than that of 7’. As a partial compensation, the denominator 
v is also greater than v’, but it may be shown that the difference in the denomi- 
nators will usually be trivial if k is small relative to n. We therefore expect r 
to be greater than 7’. Now for fixed 7; , me and significance level, it is well 
known that the power function (28) is monotone increasing with r. Hence, 
other things being equal, the increase in 7 due to the covariance variates leads to 
a more powerful test. 

The two power functions, however, differ in another respect, in that with co- 
variance the value of nz is reduced from (n — p+ 1) to(n —-p—k+ 1). This 
decrease in the number of degrees of freedom in the denominator of F will to 
some extent offset the gain from an increased 7. Examination of Tang’s tables 
[10] indicates, however, that if the degrees of freedom are substantial, this effect 
will not be important. Moreover, in most practical applications, k is likely to 
be only 1 or 2. Hence, as a first approximation the effect will be ignored, though 
to do so tends to overestimate the advantage of covariance. 

Suppose now that + = rr’, where r > 1. Since 7’ is proportional to N, the 
size of sample taken from each population, we could make 7’ = 7 by increasing 
the size of sample (when covariance is not used) from N torN. This suggests 
that the ratio 7/7’ can be used, as a first approximation, to measure the relative 
accuracy obtained with and without the use of covariance. This measure carries 
approximately the usual interpretation that the inferior method would become 
as good as the superior method if the sample size for the inferior method were 
increased by the factor r. A further refinement could be made to take account 
of the difference in the nz values. By trial and error applied to Tang’s tables, 
one could determine 7’ so that the two power functions would be as nearly coin- 
cident as possible. 

In practice, the ratio 7/7’ must be estimated from the data. From the power 
function in equation (28) it is found by integration that the mean value of y is 


(27 + p)/(ne — 2), 


so that an unbiased estimate of 7 is 


4{ (nz — 2)y — p} = | @—2) p a i} 
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This suggests that the quantity 


(nz — 2) Py. 
N> 


should be calculated both with and without covariance. The ratio of the two 
values will probably not be an unbiased estimate of 7/7’, but may be used pend- 
ing further information about its sampling distribution. This type of calcula- 
tion is made for the numerical example in section 8. 


17. The case of a row by column classification. Thus far the discussion has 
been confined to the case where there are only two populations. The technique 
may also be used when there are more than two populations. The difference 
6; between the two population means is replaced by some linear function of the 
population means. As an illustration we consider a row by column classification, 
the case that arises in the numerical example. No detailed proofs will be given, 
though it is hoped that the theory can be fairly easily developed from the mathe- 
matical model. 

A typical variate is 2;,,2, where 7 = 1, 2, --- p denotes the variate, w = 1, 
2,-:- r denotes the row and z = 1, 2, --- c denotes the column, there being 
one observation in each cell. The variates x;,z follow a multivariate normal dis- 
tribution, with covariance matrix o;;.; and means 


k 
E(Xiws) = pi + Piw + Viz + 7 Bie(Vewe adil X..), 
é=1 


where p;,, denotes the effect of the row and y;, that of the column. Without loss 
of generality we may assume that 


Lo Pie = Dy Vis = 0. 
In addition, there exists a known set of variates t, such that 


Ye=th, t=O. 


That is, the column constants have a linear regression on a set of known numbers. 
The following are the maximum likelihood estimates of the relevant constants. 


k 
be = D> Ew E™, 
n==1 
where 


Ey, = x , , . +2 lz Ly-2) | 
4in — Viwz {tre = 2 qe: t. 2 > . 
os a} 


a tien, ars 2 Diz X¢-2) 
- Dit: 


6, 
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In the notation used for numerical calculation, 


2 _ (di — Do bie d:) di 
i= —. ae eee — where d; = Zz t, Xi-s 
r>, e r>, 4 ” 


the quantity X;., being the column fofal. Finally 


rCOi;.t => Ej;.t = Ei; — > bE”. 
g 


The distributional properties are similar to those in the two-population case, 
The quantities £;;.; follow a Wishart distribution ith (re — r — 1 — k) df. 
and covariance matrix o;;. The variates d; follow a multivariate normal dis- 
tribution with means ré;2¢: and covariance 


oi;.2(r2e. + SE" d; d,) = U0%;-¢ (say). 
Consequently, 
y = 2E* did 


is distributed as vpF/(re — r — p — k) with p and (re — r — p — k) df. and 
parameter 


tr = 4(r3dt2)’do0"*5,5;/v. 


Thus in the numerical example, with r = 12,c = 4, p = 2, k = 1, this procedure 
would have given an F test of the null hypothesis 7 = 0, where F has 2 and 33 
d.f. However, the contribution from 2 degrees of freedom was deliberately 
omitted from the quantities E;; , so that F actually had 2 and 31 df. 


PART III 


18. Justification of the ‘dummy variate’ approach. It remains to show that 
the method of calculation used in the example (sections 5 and 6) is equivalent to 
that derived from theory. There arg two chief points to prove. First, that the 
M’s found from the equations 


(30) dD, Eij.¢M; = di 


are proportional to the corresponding L’s found from the equations 

(31) =. Si L; = d; 

where the suffix a denotes summation over both x; and 2; variates. 
Now, since S;; = E;; + d; d;/240, equations (31) are the same as 

(32) > E.;L; = d;(1 a 7 L; d;/240). 


Hence the L’s in (31) are proportional to the values found from the equations 
(33) D Ey Lj = di. 








ad 


re 
33 


ly 


to 
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But it is well known that if the L; are eliminated one by one from equations (33), 
we obtain 


a Ey. Lj = d; ’ 
2 
which is the same as (30). This proves the first point. 


The second point to establish is that the F test in the example is the same as 
that obtained from theory. In section 15, it was shown that 


(34) | Dd E** di dj/v 

4,2 
is distributed as pF /(n — p — k + 1). In the analysis of variance of Table 6, 
section 6, the quantity follo "ing the same distribution was 


(S. — S;) 


(35) (340 = S,)’ 


where 
Se = >, S*d;d;, S; = >) S" d& d,. 

a én 
Since equations (31) and (82) have the same solution, we must have 


S” = B7(1 — 2) L;d;/240) = E%(1 — S./240). 


Multiplying both sides by d; d; and summing over all 7, 7, we obtain 
S, = E.(1 — S,/240) = E,(1 + E,/240), 
where E, is defined analogously to S.. Similarly 
S; = E,/(1 + E;/240). 
Hence 
Soa-S; E,-EK E,— 


™ HO—S. MOF ER, 0 


Transform the variates x; , x; into variates x; , x;, where 2; = 2; — Dba; . 
It is easy to see that this transforms 
D EV did; into DOE" ded, + Do EY di d;. 
a tj 
That is, 
EB, = Ey + 2) E** di dj, 
ij 


since the quantity on the left is invariant under non-singular linear transforma- 
tions. Hence from (36), 


(Sa _— S:) ai reget / yf 
(240 — Sa) y oa 
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From (34) and (85), this establishes the equivalence of the F tests. While the 
proof has been given only for the type of data encountered in the example, the 
same method will apply to other types of data. 

In conclusion, we wish to thank the referees for many helpful suggestions in 
connection with the presentation of this paper. 
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ON THE KOLMOGOROV-SMIRNOV LIMIT THEOREMS FOR 
EMPIRICAL DISTRIBUTIONS 


By W. FELLER 


Cornell University* 


Summary. Unified and simplified derivations are given for the limiting forms 
of the difference (1) between the empirical distribution of a large sample and the 
corresponding theoretical distribution and (2) between the distributions of two 
large samples. 


1. Introduction. .Let X1,---, Xw be mutually independent random vari- 
ables with the common cumulative distribution function F(x). Let X : yet 
Xx be the same set of variables rearranged in increasing order of magnitude. 
The empirical distribution (or sum-polygon) of the sample X,,--- , Xw is the step 
function Sy(x) defined by 


0 forz < X? 


(1.1) Sy(z) = ; for X, <au< Ren 


1 forz > Xp». 


In other words, N- Sy(x) equals the number of variables X, which do not exceed 
z. We expect intuitively that Sy(r) — F(x) as N — «. In fact, if this were 
not so the notions of distribution and sample would be meaningless. The so- 
called w*-criterion of von Mises [4] provides rough estimates for the probable 
deviations of Sy(x) from F(x) for certain forms of F(x) (cf. von Mises [4]). A 
much stronger result is due to A. Kolmogorov and is of great interest in the 
theory of non-parametric estimation (Kolmogorov [3]). The maximum of the 
deviation | Sy(x) — F(x) | is a random variable Dy whose distribution is easily 
seen to be independent of the special form of F(x) provided only that F(x) is 
continuous.’ The exact distribution of Dy is not known, but Kolmogorov found 
that N*Dy has a limiting distribution. More precisely we have 

THEOREM 1 (Kolmogorov [l]). Suppose that F(x) is continuous and define 
the random variable Dy by 


(1.2) Dy = l.u.b.| Sy(z) = F(x) |. 

*Research under an ONR contract. 

1 This fact will not be used explicitly in the sequel but follows as a byproduct from our 
proofs. A simple direct proof consists in considering the random variables =, = F(X) 
which are uniformly distributed; the maximum deviation Dy of the empirical distribution 
of the new sample {=} from the uniform distribution has the same distribution as Dy; 
ef. Kolmogorov [1]. 
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Then for every fixed z > 0asN > ~ 


(1.3) Pr {Dy < zN™} — L(z) 


where L (z) is the cumulative distribution function which for z > 0 is given by either 
of the following equivalent relations” 


; . ™ 

(14) L@)=1-2>0 (-1"e"* = (2e)'s* Dee 
v=1 v= 1 

For z < 0 we have, of course, L(z) = UV. 

Equally interesting is Smirnov’s result concerning the maximum difference 
between the empirical distributions of two samples with the same cumulative 
distribution. 

THEOREM 2 (Smirnov [5]). Let (X,,--- , Xm) and (¥,,--- , Yn) be two sam- 
ples of mutually independent random variables having a common continuous dis- 
tribution F(x). Let S»(x) and T,(x) be the corresponding emyirical distribution 
functions and define a new random variable Dn» by 


(1.5) Dan = \.u.b.| Sx(2) — T(x) |. 
Put 


1 a V = ment 
(1.6) men 


and suppose thatm —» ©, -> ~ so that 


= m 
) , 


(1.7 


— a, 
n 


where ais aconstant. Then for every fixed z > V 
(1.8) Pr {Dn < 2N~*} — Le), 


where L(z) is the same as in (1.4). 

The original proofs (Kolmogorov [1] and Smirnov [6]) are very intricate and 
are based on completely different methods. IXolmogorov’s proof is based on an 
auxiliary theorem of equal depth proved in a separate paper (Kolmogorov [2)). 
An alternative proof of Kolmogorov’s theorem is due to Smirnov [5]. However, 
Smirnov derives both theorems as corollaries to much deeper (but less useful) 
results concerning the number of intersections of the graphs of Sy(x) and F(x) + 
eN and of S,,(x) and 7,,(x) +. nN, respectively. It is, therefore, not 
surprising that Smirnov’s proofs require a powerful technique and many auxiliary 
considerations. It is the purpose of the present paper to present unified proofs 
of the two theorems which are based on methods of great generality.’ The new 

2 The equivalence of the two formulas in (1.4) is a well-known relation often called trans- 
formation formula for theta-functions. We shall only prove the first representation in 
(1.4). The second is more useful for small z. A table of L(z) is given in Smirnov [6]. 
It is reprinted in the present issue of the Annals of Mathematical Statistics (pp. 279-281). 

3 Among other results which can be proved by the same method are certain limit theo- 
rems for ruin and first-passage time problems in the theory of diffusion and random walks. 
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proof is not simple but simpler than the original ones. At any rate, it requires 
essentially only routine manipulations with generating functions and their 
limiting form, the Laplace transforms. However, the paper aims mostly at a 
unification of methods. 

As a byproduct of the proof we obtain 

THEOREM 3. Let Ay be the number of points x where the step-polygon Sy(zx) 
of Theorem 1 leaves the strip F(x) + zN~*. The expected value of the random vari- 
able Aw satisfies the asymptotic relation 


(1.9) E(Ay) ~ 2(2nN)*{1 — &(22)}, 


where ®(z) ts the normalized Gaussian distribution. 

An analogous corollary to Theorem 2 was given by Smirnov [8]: formula (1.9) 
holds also for the number of intersections of the graph of S,,(x) with the step- 
polygons 7',(%) = zN~*. These results should come as a surprise to most statis- 
ticians. According to Theorem 1 there is a positive probability that Ay = 0 
and nevertheless E(Aw) is of the order of magnitude N*. The explanation 
lies in the fact that if S,(x) crosses the curve F(x) + zN ~ at some point then it 
is extremely likely that S(x) will in some neighborhood continue to fluctuate 
around values F(x) + zN~?, crossing that curve a great many times. The differ- 
ence Sy(x) — F(x) exhibits, in the limit N — ©, many small oscillations. This 
phenomenon is related to the well-known fact that the path of a particle subject 
to the Einstein-Wiener diffusion process has no derivatives. 

Instead of the absolute values of the differences we may consider the differ- 
ences themselves and derive two parallel theorems for the maximum and the 
minimum. As an example we shall prove 

THEOREM 4. With the notations and assumptions of Theorem | let 


(1.10) Dk = 1.u.b.{ Sy(x) — F(x)}. 
Then 
(1.11) Pr{ Di < zN} 31 -— &™. 


The proof is simpler than that of Theorem 1 but uses the same method. 


2. Notations and preliminary remarks. For printing convenience it is desir- 
able to avoid complicated subscripts and we shall therefore use the following 
notation for bincmial coefficients 


(2.1) C(n, k) = (”) , 


Similarly, for the general term of the binomial distribution we shall write 


(2.2) Bn, k; p) = C(n, k)p*(1 — p)"™. 
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If A is an event, A will denote its negation (complementary event). Finally 
(2.3) Pr{A | B} 


denotes the conditional probability of A for given B. 

Our proofs depend on a special case of the continuity theorem for charae. 
teristic functions. Since we shall deal only with probability density functions 
f(® which vanish for t < 0 it is preferable to use, instead of the characteristic 
function, the Laplace transform 


(2.4) si) I ” eft) alt 


(This amounts to using the variable — s instead of the usual zs and therefore 
¢(s) obeys the formal rules for characteristic functions. ) 

For any sequence {u;}(k = 1, 2, ---) of non-negative numbers we define the 
generating function u(r) by 


(2.5) u(A) = : ur. 


Now let 6 > 0 be fixed and consider the step-function f;(¢) defined by 
(2.6) felt) = ue for (k — 16<t < ké 
(kK = 1, 2,---; fs) = O fort < 0). Its Laplace transform is 


(2.7) ¢3(s) = i u(e*). 


We have, therefore, the continuity-theorem: If, as 6 — 0, 
(2.8) du(e")  ¢(s), 

then for every fixedt > 0 

(2.9) up —> f(t) when ki >t; 
conversely, if (2.9) holds then (2.8) is true. 


3. Proof of Theorem 1. Since F(z) is continuous it is possible to define num- 
bers 2, such that 


(3.1) F(t) = a ‘ciaitadil 


This definition is unique except when F(x) = k/N within an entire interval, in 
which case we define x; as the left endpoint of that interval. 

Let c > O be an integer. We shall evaluate the probability of the event 
Dy > c/N and we shall later put 


(3.2) c = zN', | 
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Suppose first that for some particular zx 
c 
(3.3) Sy(x) — F(x) > N° 


This point x is contained in a maximal interval in which (3.3) holds and at the 
right endpoint é of this interval we shall have 


(3.4) Sy(t) — F@ = x: 

Now Sw(é) is necessarily a number of the form r/N with an integer r. Since c 
is an integer also F(E) = k/N and hence ¢ = x for some k. From (3.4) we 
conclude that 


(3.5) Xize< te, Xbyey > & 


or in other words: exactly k + c among the N variables X, are smaller than 2, . 
Denote this event by A;(c). The inequality (3.3) takes place for some z if, and 
only if, at least one among the events A,(c), --+ , Aw(c) occurs. The argument 
applies equally to c < 0 and shows that the event Dy > c/N occurs if, and only tf, 
at least one among the events 


(3.6) A,(c), Ai(—e), As(c), Ax(—c),--- , Aw(c), Aw(—c) 


occurs. 

Let U, and V, be the events that in the sequence (3.6) the first event to occur 
are A,(c) or A,(—c), respectively. More formally, the events U, and V, are 
defined by 


r= A;(c)Ai(—c) eee A,-1(c)Ar-1(—c)A,(c) 
V, = Av(c)Ai(—c) - ++ Ay-1(c)Ay-1(—c)A,-(c) Ar(—c). 


These events are mutually exclusive and therefore 


(3.7) 


(3.8) Pr { Dy > ct = >) Pr {U,} + a Pr {V,}. 


N 


From the very definitions we have the following two fundamental relations 


r=1 


k 
Pr {Az(c)} = 2» Pr {U,} Pr {A.(c) | A,(c)} 
+ 2» Pr {V,} Pr {Ax(c) | A-(—c)} 
(3.9) ' 7 
Pr {Ai(—c)} = dX Pr {U,} Pr {Ax(—c) | Ar(c)} 


+ dX Pr {V,} Pr {Ax(—c) | A-(—o)}. 
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This is a system of 2N linear equations for the 2N unknowns Pr {U,} and 
Pr {V,} and we proceed to solve it by the method of generating functions. 

By definition of x, we have Pr {X, < a} = k/N. The probability of the event 
A;(c) (that the same inequality holds for exactly k + c different v’s) is therefore 
given by 


(3.10) Pr {A,(c)} = BIN, k + ¢;k/N) 
(ef. (2.2)). Similarly, it is readily verified that for r < k 


(3.11) Pr {Ax(c) | A-(c)} = BIN —r—c,k —717;(k — r)/(N —1)). 

and 

(3.12) Pr {A.(c)|A-(— c)} = BIN —r+c,k — r+ 2c; (k — r)/(N —1)). 
The last three equations hold also fore <0. They can be written in a more con- 
venient form in terms of the quantities 

pkte 


(3.13) + pile) =e 


In fact 


(3.14) Pr {Az(c)} “ pi(c) pri —c) 
p(0) 
(3.15) Pr {A;(c) | A-(c)} = pr—r(0)pr—x(—e) 


Pr-r(—c) 
9) \ a . 
(3.16) Pr {Az(c) A,(—c)} = pr—r(2c)pw—e( Cc) 
Px-r(C) 
If these expressions are introduced into (3.9) the second factor in the numerator 
cancels. A further simplification is achieved on introducing new sets of un- 
knowns 


pn(0) 
Pw-r(—C) 
The fundamental equations (3.6) then reduce to 


k 


k 
pile) = Dy Ur pr-r(O) + Do 0, p-r(2c) 
1 r==1 


3.17) - = P U; ae : 
(3.17) u r {U,} a 


(3.18) 


& k 
pi(—c) = dX Ur Di-r(— 2c) + 2d Vy Dr-r(O). 


This system is of the convolution type and can therefore be solved by means 
of generating functions. The essential point is that the p,(c) are defined for 
all k and that the system (3.18) therefore determines the unknowns u, and » 
for allr > 0. We put 


(3.19) ati a ‘ ai 0m : ve dt 
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id ill 
. (3.20) p30) = N*Y plo) a. 
re k=1 


(The factor N™ + serves to simplify formulas.) Then obviously 
P(A; c) = u(A)p; 0) + v(A)pQs 2c); 





(3.21) 





p(aA; —c) = u(r)p(a; —2c) + v(A)p(; 0). 


From here we find u(A) and (A). Equation (3.17) then determines Pr { U,} and 
Pr {V,}. Actually we are interested only in the two sums occurring in (3.3). 
We put 







(3.22) & = =o 2 > Pi-r(— C)Ur ’ mm > mm 2 > Pr—r(C)U, . 


Again, the & and 7 are defined for all k (alsok > N). From (8.17) we have 










(3.23) dX Pr {U,} = &y, a Pr {V,} = ow, 








and hence finally, by (3.8) 
(3.24) Pr {Dy > c/N} = ty + ow. 


In (3.22) we find again simple convolutions leading to products of the corre- 
sponding generating functions. Thus 






. eat = u(r)p(a; —c)N} 














- &(A) = (0 
‘oo : (a) “ < 
i k _ V\A)P\A; C 
n() = dX mA = - ma(0) ; 





We now pass to a study of the limiting form of these generating functions 
as N —> o andc— o in accordance with (3.2). Consider a fixed ¢ > 0 and sup- 
pose that 








(3.26) 5 sabi 





From well-known properties of the Poisson distribution it follows then that 
(3.27) N*px(c) — (2at)exp(— 2°/2%). 


Accordingly, the continuity theorem of section 2 implies (as can be verified dir- 
ectly) that 






(3.28) Pe "5 2N*) — (22) | v7 exp (—ts — 2*/2t) dt 
? 0 





= (2s) exp(—(2 322)4), 
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(the last integral is well known and can be evaluated by elementary methods, 





the square-root is always positive). We see in particular that the limiting form ( 
is the same for p(\; c) and p(A; — c). It follows therefore from (3.21) directly ; 
that ! 













9.-2\3 
39 , -IN) — tie a(n) — __ exp(—(2sz")') 
mae = mn . wn 1 + exp (—(8s2?)!) 


Using this and the fact that py(0) > (2nN)~* we conclude from (3.25) that 


lim N7e(e""*) = lim J 7 ne) 


N-o2 





(3.30) _ (@" _ exp (- (8sz*)") — 
2s 1 + exp (— (8sz?)!/?) a 
Expanding ¢(s) into a geometric series we get 


1/2 @ 
(3.31) ¢(s) = (37) a (—1)"” exp (— (8sv*2*)”?). 


From the evaluation of the integral in (3.28) we conclude that ¢(s) is the Laplace 
transform of 


(3.32) fi) = : (1) exp (—2/72/2). 


The continuity theorem of section 2 in conjunction with (3.30) and (3.26) shows 
that 


(3.33) lim ~y = lim qv = f(1). 
N-©o 


N-o 


In view of (3.24) this accomplishes the proof. 


4. Proof of Theorem 4. This proof is simpler than the preceding one inas- 
much as we are now interested only in the events A;(c) fore > 0. This time we 
define U, as the event that k is the smallest subscript for which A;(c) occurs, that 
is, U, = Ai(c)A2(c) --- A,«a(c)A,(c); no analogue to the event V, will be used. 
With the same notations as before (3.9) is replaced by 









k 
(4.1) Pr {Ax(c)} = dX Pr {U,} Pr {Ax(c) | A-(c)}, 
and hence (3.21) by 
(4.2) p(r;c) = u(r)pQ; 0). 
Here p(A; c) is the same as before, so that (cf. (3.29)) 
(4.3) lim u(e“"*) = exp (— (2s2*)"”). 


N-0 







Again, the first equation (3.25) holds without change and therefore we get in- 
stead of (3.30) 






‘ 1/2 
(4.4) lim N"g(e"™) = (3) exp (— (8sz*)"). 
N-o 
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From (3.28) this is the Laplace transform of 
(4.5) fi) = €? exp (— 22/1). 

As before we conclude that £y — f(1), which accomplishes the proof. 

5. Proof of Theorem 3. We have seen in section 3 that the intervals in which 
(3.3) holds are in a one-to-one correspondence with the events A;(c). Hence 
(5.1) E(Aw) = = Pr {A,(c)} +2 Pr {Ai(—o)}. 

To evaluate the sums we use (3.10). If N — o and again c = 2N?, k/N —t, 


then by the central limit theorem 


+ be B/W) —» oP (—# /2l — 0) 
(52) BIN, k + 6; k/N) > 


It follows then from (3.10) that 
1 e 

(53) N~"?= Pr {Ax(c)} — (22)? | {t(1 — i)}7*? exp (—2?/2u(1 — 2) dt. 
0 


Call the right hand member R. After the substitution ¢ = sin’ (¢/2) we find 


dk a —1/2 ,, —2 Patan 
_* —8(2r) “2 sin ¢ exp (—2z:°/sin’ ¢) dd 
0 


3/2 
(5.4) = 8(27)~"* z exp (—2:”) | exp (—2:’ cot? ¢) d (cot ¢) 
0 


= —2 exp (—2:”). 


Since R — 0 as z — © we conclude that 
(5.5) R=2 H exp (—2z”) dx = {1 — &(2:)}(2x)™*. 


The same asymptotic estimate holds for the other sum in (5.1), and hence Theo- 
rem 3 is proved. 


6. Proof of Theorem 2. Reorder the two samples in ascending order 
of magnitude and denote the rearranged samples by (Xi,°°+, Xn) and 
7; ,°++, ¥2). When speaking of the graphs of the empirical distributions 
Sn(z) and T,(x) we shall suppose that they have been completed by adding 
vertical segments so that the graphs become step-polygons. We shall put 


m 
(6.1) ats > 


Then, according to (1.6) and (1.7) 


(6.2) _ a, 
q 
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Without loss of generality we shall suppose that 
(6.3) p <q. 


In order to carry over the proof of Theorem 1 it is necessary to define the 
events A;(c) in a judicious manner. For every integer k > 0 let »% be the num- 
ber of variables X, which are smaller than Y;,. In other words, »% is defined 
as the integer for which 


(6.4) x < ¢ Xow 
Finally put 


(6.5) a, = | | = |? | 
n q 


where, as usual, [z] denotes the greatest integer contained in z. 
For 0 < k < n let A;(c) be the event that 


(6.6) Ve = Ak+e- 
The possibility of applying the proof of section 1 depends on the following 
Lemma. Whenever 


(6.7) a. > < >0 


then at least one among the events A,(c), Ai(—c),---, An(c), An(—c) occurs. 
Conversely, if one of these events occurs then 


an 
(6.8) Rua > (« 2) /n 


Proor. If (6.7) holds then either for some zo 
(6.9) Sm(t) — Tn(ao) > 


or the reversed inequality holds with c replaced by — c. It suffices to consider 
the case (6.9). For sufficiently large x we have S,,(z) = T,(x) = 1. Hence the 
graphs of S,,(x) and 7,(x) + c/n must intersect at an abscissa § > 2. The 
point of intersection lies necessarily on a horizontal segment of the graph of 
Sm(x) and a vertical segment of 7,(x) + c/n. Hence there exists a k such that 
¢ = Y,” and, moreover, 


(6.10) Ta —) + < < Sul) < TalE +) + =. 
This amounts to saying that 
(6.11) kK-1lte wm -k+e 

n m n 


In view of (6.3) and (6.5) this relation implies (6.6). 
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Conversely, suppose that the event A;(c) occurs and let c > 0. Put again 
* o,e 
¢= Y;,. Then, by definition, 


(6.12) S,.(g) = Sete T(t) = 
m m 


It follows that 


. 
n- 


(6.13) Sn(&) > - 


which in turn implies (6.8). This proves the lemma. 

Theorem 2 is concerned with values of ¢ such that en™' = zN~*; in passing to 

the limit we must therefore put 
(6.14) c = 2(n/p)'. 
Accordingly, the relations (6.7) and (6.8) are asymptotically equivalent and our 
lemma shows that, asymptotically, the probability of (6.7) is the same as the 
probability that at least one among the events Ai(c),--- , Aw(—c) occurs. 
To evaluate this probability we proceed exactly as in section 3. The events 
U, and V, defined by (3.7) and the fundamental relations (3.9) hold again. 
However (3.10) — (3.12) have to be replaced by new evaluations. 

It is easily seen that the probability that exactly r among the X, are smaller 
than Y; is the same as the probability to extract exactly r white balls before the 
k-th black ball from an urn containing m white and n black balls (assuming 
that all orders are equally likely and that balls are not replaced). In this way 
one finds 


(6.15) Pr {Az(c)} = 


Pr {Ax(c) | A,-(c)} 


(6.16) pom C (Gite ba Gr+e + kt em 1k oe 1)C(m +2— Ak+e — k,n a k) 
y C(m + 1 — Grae — 7,2 — 7) 


Clante +h —1,k —1)C(m+n — ays. — k, n — k) 
C(m + n,n) 








Pr {A;(c) | A-(—e)} 


(6.17) — Cee — Are HK -— 1 — lk—r—1)C(m+n— ay. — k, n—k) 
C(m + n — a,-. — 7,0 — 1) 





The second binomial coefficient in the numerator is common to the three ex- 
pressions and cancels when the expressions are introduced into (3.9). These 
fundamental relations assume a more natural form if the occurring binomial 
coefficients are enlarged to terms of a binomial distribution. It is easily veri- 
fied that the first of the equations (3.9) reduces to 


B(aiz4e +k — 1,k — 139) 
nse 
_ BlQkse — Ore HK —-—1r—1k—r—1;¢q) 
ay = iru ee ae 
Bite ~ Gre HK -—7—- 1k -—7— 15) 
Bim + 0 — Gr—e — T, N — T;Q) 


k 
+2) Pr {V;} 
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The second equation is obtained on replacing the combination k + ¢ by 
k-— ec. 
Instead of (3.17) we put 


Bim + n,n; q) 

y » Boe Re ’ ? 7 

Pr (U;} Bim +a - Grie — TN — 7; q) 
Bm+n,n3g) 

Bim + 2 — Gp — 7,2 — 73q) 


Ur 





(6.19) 





v, = Pr {V,} 


Then (6.18) becomes 


B(Qite +k — 1,k — 139) 
k 

= >) tte Batre — Orre +k —r-lk-r-—1;q) 
r=1 


k 
a Z. 0,B(dkre — Gee Hk —r—-lk—-r— 1; 9). 
r=1 


This corresponds to the first equation in (3.18). Unfortunately (6.20) is not 
of the pure convolution type since ax4- — r+. and ay4- — a,—- are not, functions 
of the two variables k — r and c. The trouble comes from the fact that a, 
as defined by (6.5), is not a linear function of k. It is, however, plausible that 
we shall commit only an asymptotically negligible error if we omit the brackets 
in (6.5), that is, if we replace a, by pk/g. Purely formally (6.20) then reduces 
to the first equation in (3.18) with 


(6.21) pi(c) = B (te —1,k-1; ‘). 


(Here the first argument in the right hand member is no longer necessarily an 
integer, and the factorials in the definition (2.2) should be interpreted by means 
of the gamma function.) To the new system (3.18) the considerations of section 
3 apply almost word for word: the only difference lies in the new norming (6.14) 
(which replaces (3.2)) and that instead of (3.26) we shall naturally let k/n > t. 
Thus the limiting form of Theorem 1 applies to the new system (3.18) with p:(c) 
defined by (6.2). 

It remains to prove that the formal replacement of (6.20) and the corre- 
sponding equation for — c, by (3.18) was legitimate. Now all coefficients 
in (6.20) are of the form B(v, r; q), and we have only changed the first argument, 
y, adding a variable quantity which in no case exceeds one unit. In passing to 
the limit we put k ~ tn and c ~ zn’p™?. It follows that we actually use only 
coefficients B(v,r;q) where v > ©,r — © and v/r—>q. Accordingly, for | 3| 
< 1 we have B(vy + 8, r;q) ~ B(v,7;q), and it is rather obvious that our system 
(6.20) is asymptotically equivalent to (3.18). 
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APPLICATION OF RECURRENT SERIES IN RENEWAL THEORY 


By Aurrep J. LorKa 


Metropolitan Life Insurance Company 


Summary. The application of integral equations to renewal theory in popv- 
lation analysis and problems of industrial replacement is beset with certain diff- 
culties which have been particularly discussed by W. Feller (these Annals 194} 
vol. 12 pp. 243-267). Some of these difficulties are avoided if the data of the 
problem are introduced into the analysis directly in the discontinuous form 
(tabulated by class intervals) in which they are usually supplied in a concrete 
case. A numerical example based on population statistics is presented, illustrat- 
ing how, using discontinuous data, a recurrent series takes the place of the integ- 
ral equation, and a finite exponential series appears in place of the Heaviside 
expansion of the previous solution. There is close analogy with the procedure 
previously presented, but with factorial moments appearing in place of ordinary 
moments. 

The fundamental data being given for values of the replacement function at 
discrete intervals only, some question arises as to the applicability of the solution 
as an “interpolation” formula for non-integral values of the time /, and as to the 
effects of subdividing the class interval of the original data. 

In the actual computation of the factorial moments a shift of origin by one- 
half class interval becomes necessary. An algorithm for effecting this shift is 
presented. 




































1. Methodology: Alternatives Available. 


All application of mathematics to concrete situations involves a greater or less 
degree of conventionalisation, a substitution, “in place of intractable reality, of 
an ideal upon which it is possible to operate.’” 

This conventionalisation may be only such as to do little violence to the con- 
crete data, as for example when, dealing with a large population, we treat the 
number N(/) of individuals at time ¢ as a continuous variable, knowing perfectly 
well that strictly speaking it varies by jumps of one unit at a time.” 

In dealing with any particular concrete case there may be considerable choice 
as to the mode in which the conventionalisation or idealisation is carried out, 
and the particular place or step in the scheme at which it is introduced. A good 
illustration of this is met in the treatment of renewal theory, as applied to human 
populations or other biological or industrial aggregates. 

The majority of authors who have dealt with the subject have set up their 
aenieenane oquiime in terms of continuous variables. Many have gone fur- 


1 Nature, Vol. 110 (1922), p. 764. 
2 If the susidhilien is subject to extreme variation in numbers, such that N(t) passes 
through small values, this disregard of their discontinuity may not be permissible. 
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ther than this in the process of conventionalisation and have assumed for the 
renewal function (net reproductivity) some more or less appropriate mathemat- 
ical expression, such as a Charlier or a Pearson [1] frequency distribution, and 
have, wherever possible, carried out by standard methods the integrations in- 
volved. 

Others, while retaining the formulation of the fundamental equations in con- 
tinuous (infinitesimal) form, have made no specific assumptions regarding the 
analytical form of the renewal function, and have carried out the numerical in- 
tegration by one of the established methods available for the approximate in- 
tegration of arbitrary functions. 

But there has also been a minority of authors who deemed it most appropriate, 
since the data of the problem are actually furnished in tabular (and hence dis- 
continuous) form, to apply from the start discontinuous methods in formulating 
the fundamental equation for the problem. This equation then defines a recur- 
rent series. 

The most recent and also the most concise exposition of this approach to the 
problem is a paper by W. Dobbernack and G. Tietz presented at the Twelfth 
International Congress of Actuaries, 1940, Proceedings, vol. 4, p. 233. These 
authors, however, do not give any numerical application, and in consequence 
certain aspects of the analysis are not touched upon by them. A more detailed 
presentation, including numerical applications, was given by the late S. D. Wick- 
sell’ who, however, used only roughly approximate data (an over-all average net 
reproductivity for ages 20 to 44) and also introduced certain linear interpolations 
which would not be appropriate with more exact data, and which become un- 
necessary in the numerical operations if moments are introduced as indicated in 
what follows. “ 

The purpose of the present paper is to exhibit this modification of the method 
of recurrent series, and at the same time to illustrate its relation to the method 
which proceeds in terms of a continuous variable, leading to an integral equation. 

The B(t — a) women born in the calendar year (¢ — a), that is, between the 
times (¢ — 3 — a) and (¢ + 3 — a), will be a years old some time during the 
calendar year /, that is, between t — 4 and ¢ + 3. If their births were evenly 
distributed over the year ¢ — a, so will their birthdays of age a be over the year 
t, and their average age during that year will be a and the average number of sur- 
vivors to that age during the year ¢ will be approximately B(t — a)p(a), where 
p(a) is the probability, at birth, of surviving to age a. If the annual female 
reproductive rate, counting daughters only, is m(a) at age a, then the B(t — a)- 
p(a) survivors will, during the calendar year t, give birth to B(t — a)p(a)m(a) 
daughters. If B(é) is the total number of births of daughters in the calendar 
year /, then evidently, for positive values of 1, 


(1) Bit) = > B(t — a)p(a)m(a), 


3 [2]; see also (3). 
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or, to simplify the notation, 


@ 


(2) Bit) = D> ca Blt — a). 


Equation (1) or (2) defines a recurrent series of the general form 
(3) B(t) = Bt — 1) + @Bt — 2) + +++ + Bt — w), 


where some of the coefficients c may be zero and where w denotes the upper limit 
of the reproductive period. 
The trial substitution 


(4) Bit) = Qu" 
in (3) gives 
(5) l=oar+ or +---+ cer”. 
The substitution (4) therefore satisfies (3) provided that z is a root of the equa- 


tion (5) of degree w for x; and the same is evidently true for the more general 
substitution 


j=w 


(6) Bit) = » Q3x;', 


where x; , with 7 = 1, 2, --- w, are the w roots of (5). 

Equation (5) leaves the w coefficients Q; indeterminate. In general they ap- 
pear as arbitrary constants. In any concrete application they may be deter- 
mined by “‘initial’’ conditions; that is, in order to make the problem determinate, 
it is necessary to be given the values of B(é) for w successive integral values of t, 
or some equivalent data. 

While, for convenience in description, the analysis has been developed in 
terms of the year as time unit, the formulae are evidently independent of this 
choice of unit, provided that the unit employed is adequate for practical appli- 
cation. 

Whatever the unit employed, for the direct application of (1) and (3) to a con- 
crete case it is necessary to have the data in such form that values of p(a)m(a) 
are known for integral values of a. The pertinent statistics do not usually come 
in that form, the fertility being usually known only for five year age groups, and 
though it may be sufficient for practical purposes to regard these quinquennial 
values as representing p(a)m(a) for the midpoint of the group, this yields p(a)- 
m(a) for fractional values of a, as measured in five year units. We may then 
proceed as follows: putting 


(7) r=l+y 
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in (5) this becomes 
={atateot--: +0¢,} 
+ [ar + 2c, + 33 + --- + wea}y 
+a + 8a + 6a ++ +82 OD Se 
+ {ot des + 105, +--+ +82 TOW?) hy 
ie ein 
+ fes}y” 
-> > r 7 Chey 


In application to a particular population, we shall usually have the condition 
=Q for a=1,2,-°-<e 
where a is the lower limit of the reproductive period. 
The expressions in brackets (coefficients of successive powers of y) will be recog- 
nized as cumulations S, of the values of the function c, , summed backwards to 


the “diagonal” element c; , where h is the exponent of y. In terms of moments 
m of the function ca , equation (8) can be written 


m3 — 3m, + 2m 
or, using the symbol m,,; to denote the hth factorial moment, equation (9) takes 
the simple form 


— —t 


(9) = m + my Fi al. ye + yt ees bay” 


(10) T a y’ 


h=0 
In these expressions the moments m, a mx are those taken about a = 0. 
Actually, the net reproduction rates are given for “‘semi-values” of a, that is, for 
values of a which are odd multiples of $ (using five years as the time unit). By 
cumulation of these given values namie m, and mt, about a = —} are 
obtained.* From the latter the corresponding functions of the moments about 


a = 0 are obtained by the transformation formulae’ 


k=h (— 4)" 


/ 
Mh] a Mth—k] 


ailUKSS CUCU — 
_ 1) [kl 
Ss, = ws ( = x. i 


k=0 


(11) 


‘In these cumulations zero values of ca for 0 < a < a must not be omitted. 

5In accordance with a customary notation the symbol (—3)] denotes the continued 
product — 3}(—3 — 1)(—3-— 2)... (-3 —k +1). Inthe computation of successive terms, 
in the sums in the right-hand member of (11), by appropriately laying out the work, ad- 
vantage is taken of the fact that values of (—})[*l/k! for k = 2,3...are obtained each 
from the preceding by multiplying successively by 3, 3, 3, etc., and taking care of the sign, 
80 that fractions with complicated numerators and denominators are avoided. 
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It will be recalled that in the treatment of the problem of replacement by 
means of an integral equation,’ a solution in the form 


(12) Bit) = 2Qix;* = BQje", 

is obtained, in which the exponential coefficients r; are the roots of the equation 
(13) 1= [ e * p(a)m(a) da = [ x° p(a)m(a) da, 

i.e. 

(14) L=m—mrt+ orate =D (- 1) > 


in close analogy to equation (10) for y, with the distinction however that in (10) 
the factorial moments take the place of the ordinary moments of (14), and that 
the series in (10) is finite, terminating at the term in y”. There is also an impor- 
tant difference between the characteristic equation (13) and its analogue (5), 
namely that (5) may admit of negative roots for x, whereas (13) does not admit 
negative values for z. 


2. The constants Q. These are determined by initial conditions, as follows. 
Equation (2) can be written 


a=W a=t—1 
Bit) = Do caB(t — a) + a caB(t — a) 
a=t a= 
(15) a=t—1 
= F(t) + Zz. ca Blt ae a), 
a=1 
with 
F(t) = D> ca Blt — a) 0<iti<ae 
(16) and —_ 
F(t) = 0 t>w 
The values of B(t) being given for integral values of t, from t = —(w— 1) to 


t = 0, it can be shown that’ 


Mi 


; Feo 
(17) Z°- —— 


gic 
io 
& -_ 


Ma X5 


8 
ll 
= 


6 For a discussion of the limits of applicability of this method See [4]. 
7 The reasoning is essentially the same as in the treatment of the problem by integral 
equations. See [5] and [2, p. 39 et seq.]. 
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by In the special case that we are tracing the progeny of an initial population all 
born at the same time, say B(0) births occurring at ¢ = 0, so that 
(18) B(—1) = B(—2) = --- = B(—fw — 1) = 0 

hn the expression for Q;, in view of (5), reduces to a particularly simple form. 


For if we write the summation in equation (16) in expanded form, we have 
Fd) = 
c,B(0) + cB(—1) + c3B(—2) + cB(—3) + --: + cB(—w — 1)’ 





FQ) = o2B(0) + cB(—1) + cB(-2) + ++ + ¢B(-a — 2), 
(19) F(3) = c.B(0) + «B(—1) + --- + ¢.B(—w — 3) 
10) . ; ; 
hat 
or- 
5), sd 
nit F@) = CB (0). 
If now B(—1), --- , B(— w — 1), all vanish, then 
ws (20) X Fs! = BO){ax + a2 + ea" +++ + 602") 
(21) = B(O) by (5). 
Hence, 
B(0) 
Q; = a=w . 
(22) sual 
In particular 
Sw _e., Sl 
(23) B(O) = 2d, Q; = B(O) 2, Ten 
so that 
to S| 





(24) X a 1. 


The constant B(0) here evidently functions essentially as an arbitrary unit 
of annual births, and may with this understanding simply be put = 1, thereby 
simplifying the notation. This has been done in what follows, where con- 
venient, especially in the table of constants, Table 3 of the numerical illustration. 

The denominator in (17) or (23) can be evaluated for any root x; of (5) by direct 
summation if the coefficients c, are given or have been computed (as indicated 

ral below) for integral values of a; or, in a manner similar to that employed in passing 
from equation (5) to (8), the denominator can be expressed in terms of the cor- 
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responding roots y; = z; — 1 of (8) or (10), the cumulations of c, being replaced 
by cumulations of ac,. With the denominator so expressed, the constants 
Q; take the form, in obvious analogy to equation (9): 



























t=w 


& x’ F(t) 
(25) Q; = —— = tating 
. m3- ~~, M4 © = — 3m; + Qms y; 


m, + my; + ——— yj, + - 7 yy tt: 


2! 














The alternative procedure, to which reference was made in the preceding para- 
graph, is to operate upon the moments m,,) (taken about the origin O) by a 
process the inverse of cumulation—which we might term decumulation—and 
in this way to obtain from them the coefficients c,. The polynomial Zac,r} 
can then be evaluated directly. 

The decumulation is readily carried out by an algorithm which suggests itself 
from the schedule of cumulation. Analytically the relation between the two 
processes is expressed by the reciprocal sets of transformation formulae: 


Cumulation 
k=o—h * 
(2) Tr 3 (Ca Jawa Ss 
Decumulation 
k=w—h 
a __1)* h + k Mih+k) 


3. Constants (2 associated with complex roots x = ¢ “**’. 


The complex roots z; give rise to oscillatory terms which, in the special case 
of the progeny of a cohort of B(0) births, take the form® 


2B(0)e“* 


(28) G2. 4. H?2 


{G cos vt — H sin vt}, 


where 


8 
t 
£ 


(29) G= ac,e “* cos va 


® 


and 






a=W 


(30) H = > aye“ sin va. 


a=1 





These constants may be evaluated directly in this form, or, putting y = & + 7 
in the denominator of equation (25), they can be expressed in terms of the 
roots y; and the factorial moments obtained by cumulation of ac, .” 








8 The development of these formulae is analogous to that followed in the treatment of 
the problem by integral equations. ‘See [6]; for the more general case see also [7]. 

® The procedure in this case will be analogous to that followed in the development of 
equations (90) and (91) in [6]. 













APPLICATION OF RECURRENT SERIES IN RENEWAL THEORY 197 


(a) Numerical Illustration. For convenience and to furnish the opportunity 
for comparison, the same data (United States 1920) were here employed as in 
the writer’s earlier publications in which the problem was treated by the appli- 
cation of an integral equation. 


(b) Cumulation for values of m,. The two operations, of (1) cumulating the 
values of c, given for semi-values of a; and (2) allowing in the cumulated results 
NET FERTILITY 

m(a) 


pla) 
35 


p()= probability for a newly -born female to 
attain age ‘a’ 

m(a)=probable number of female births to 
a female between ages ‘a-+' and ‘a+4" 
(age in five year units) 


t 
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AGE IN 5 YEAR UNIT 


Fig. 1. Net Fertility p(a)m(a) White females, United States, 1920 


The verticals drawn in full and centered at mid-ages represent the original data; those 
drawn in dashed lines and centered at integral ages are interpolated. 


for.a shift of origin froma = —4toa = 0, can be conducted in one schedule as 
in Table 1. Cumulation is first carried out in the usual manner from the bottom 
line to the diagonal, with the result appearing immediately below the diagor al. 


From here on the procedure is as in the following example: Starting at the lower 
right hand corner, we find 


00780 X (—%) = —.00390 


12770 X (—}4) = —.06385 —.06385 x (—3) = .04789 
97395 X (—}) = —.48698 — .48698 X (—?) = 36254 


36254 X (—%) = —.30437. 

















a in 5-year 
units 
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a in 5-year | 


| 


| —12. 


.00000, 
.00000 
.00040) 
09630) 
31255) 
31025 
.23170) 
15090 
.05795 
00615 





m {6} /6! 


(9) 


26311 
.77790 
.41964 
—9.97080 
12.61050 
20548 
12. 38595 
3.45870 
.74135 
10520 
.00720 
.00015 


@ | 


.00015| 


.72500 | 


m {0} 


1.16635} 6.64127} 1664550} 24.34106 23. 16864 


1 
1 
1 
1 
1 


(3) 


16635 
. 16635 
16635 


16595) 
| 
06965) 


75710) 
44685 
21515) 
.06425| 
.00630 


.00015 


m7 


(10) 


1 
1 
—<4. 
8. 
—10. 
9. 
—§. 
4. 


)/7! 


.99717 
. 24432 
. 62974 


7768 
72445 
50875 
15411 
19298 
31260 


. 85390 
. 11255 
.00735 
.00015 - 


mii} 


a) | 


58318 
22445] 


72540, 
55945) 
48980) 
73270) 
28585 
07070! 
00645 
00015) 


Ke nwowr on 


at 
m {8} /d: 


an) 


. 36404 
. 22905 
. 51333 
.47121 
—7.85201 
9.19516 
—7.62843 


4.64474 
— 2.15630 
.97395 

. 12005 
.00750 
.00015 


—l 


+ 
~ 
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TABLE 1 


for integral values of age a.* 


L 


(5) 


m (2) /2! 


43738 
—3.61223| 


8.87050| 18. 14430 ~ 33.62800! 


5.14510. 
2.58565 
1.09585 
36315 
07730. 
00660: 


m{9}/9! 


(12) 


03483 
21633 
1.41875 
—4,15184 
7.19768 
—8.27564 
6.67488 
—3.87062 
1.61723 
48698 
12770 
00765 
00015 


.00015) 


} 
m(3}/3! 
' 


(6) 


— .36448 


| 
2.70917; 


05810} 19.82035 —9.91018 
.89175) 13.76225) 31.90655) —15.95328 


9.27380 
4.12870, 
1.54305) 
.44720 
.08405 
.00675)| 
.00015) 


| 
mio} /10! | 


(13) 


. 20551 

— 1.33993 
3.89235 
— 6.68356 
7.58600 
— 6.00739 


3.38679 
—1,34769 
36524 
— .06385 
.00780 
.00015 


.00127 


m4} /4! 


(7) 


.31892| 
— 2.25764) 
7.43264) 


15.48370! 
6.20990) 
2.08120! 

53815) 
.09095 
.00690) 
00015! 


| 
a 
} 
| 
| 





m s io ° 
Computation schedule for values of — = S; of net productivity function p(a)m(a) = ¢, 
J: 


m {5} /5! 


15.05650 
— .28703 
1.97544 

— 6.19387 
11.96496 
—16.81400 


24.41095 
8.92725 
2.71735 

. 63615 
-09800 
00705 
.00015 


m{unj/11! | 

| 
00001) 
19617 
27293 
67611 
26584 
04414 
50677 
—3.04811) 

1.17923 
. 30437 
.04789 
.00390 


Rs 
=F 
6. 
=—t. 
5. 


ee | 
3: 


00015 


Factor 





—21/22 
—19/20 
—17/18 
—15/16 
—13/14 
—11/12 
—9/10 
-7/8 
—5/6 
—3/4 
-1/2 


* Figures immediately below the diagonal, obtained by cumulation from the bottom 
upward of the data in Column 2, are factorial moments abouta 


Figures in the top 
line are factorial moments about a 


0. For use of factors in the last column see text. 


The several columns are thus completed, and by addition, in each column, of 
the item immediately below the diagonal, and of all the items above the diag- 
onal, the figures in the top line are obtained. These are the coefficients of equa- 
tion (10) for y. 
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(c) Decumulation. While it is not necessary to carry out the decumulation, since 
the entire computation can, if desired, be carried out in terms of y’s and m’s, 
there is a considerable interest in noting the values c, for integral values of a 
which result from the decumulations of the m’s. These, together with the 
original values for semi-values of a, are shown in Table 2 and Fig. 1. 


TABLE 2 
Values of ca = p(a)m(a) 
(1) for semt-values of a; original data. 


(2) for integral values of a; computed by cumulation of original data, shift of 
origin, and decumulation. 
a 5-year units Ce a 5-year units! a 5-year units! Ca 


0 ’ , .10607 
| 0 05795 
| O* .02268 
| 0 .00615 
0* 00116 

.00040 00015 


i | .02073* .00001 
3.5 .09630 


SIND HD Or or eS 
aoonono uo 


oo eas I soil eects oe 


*The value of c. came out negative, namely —.00570, and the value of c: came out 
+.00014. In the computation of Zac,xt these two values were arbitrarily adjusted to 
11 


zero, and c; was diminished from .02118 to .02073 to make the total 2 c. = 1.16635, sum- 


t=1 


ing only for integral values of a. 


4. The roots of equations (5) and (8). 


From the prior study already cited, the real positive and three pairs of complex 
roots for r of the characteristic equation 


(31) [ x p(a)m(a) da = [ e  p(a)m(a) da = 1 


were known. These were used to indicate the approximate location of the roots 
of (5) or (8), and more exact values were then obtained by Newton’s method of 
successive approximation. Table 3 shows the values of wu, v, etc., corresponding 
to the new roots 


z-l 
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obtained through equations (8) or (10); and, for comparison the corresponding 
values obtained in the previous publication from equation (13). The same 
table also exhibits the remaining roots and values of the constants Q, G, H. 





TABLE 3 
Constants of the series solution (6) of equation (8), corresponding to the five real and 
three pairs of complex roots of the characteristic equation (5) 

(United States, white females, 1920) 


Constants‘! | Five Real Roots Three Pairs of Complex Roots 











A. Computed on basis of recurrent series 








% | oamis*|—1.7644 | -3.812¢ |-17 17- -4, .3t/- —.19800| — = 44729) =| 47587 
° 0. 0. 0. . | 1.06498] 1.57000} 2.404 
G b pane 7.73354|—1255.04 | o | @ | 5 28003, 10.45809| 7.73103 
H 0 | 0. | 0 | 0 | 3.03939|—3.66726) 2.00874 
G/(G2+H?) : mi aoa 00080 ® | @ | 14241)  .08515| «12117 
_AKG+H) (0. | 0. | 0 | 0 | .08177| —.02986 .03148 


— 





B. Computed o on n basis of ‘integral equationt 





1 | | 
u | 02714 | | |— .1930 30 | - — .43655'— .4902 
v 0. | | 1.0724 | 1.5771 | 2.44245 
G 5.64514 | 5.15351) 10.22495, 7.40154 
H 0. 2. > 98757 — —3.72741| 3.45312 
: G/(G2+H*) | .17715 14525) —.08620; 11095 


ew 0. .08420. —.03135 = .05175 








@ tin on year emilee 

) Not computed 

*u, = log. to = —log, .97322 = .02714 
¢ Values of zx 
t See [6, p. 899] 





To determine the remaining four roots, the product of the factors (y — 4) 
(y — ye) ++: (y — yz) was divided out of the polynomial of equation (10), re- 
jecting the remainder and leaving a fourth degree equation 


y* + 120 y* + 2590 y® + 14617 y + 23118 = 0 


In the subsequent work it turned out that the roots of this were all real, and 
they were computed by obvious methods. Their values are also shown in 
Table 3. For the two numerically largest roots great accuracy was not at- 
tempted. They introduce terms with very rapid damping and presumably 
very small values of Q.” 





10 The divergence is due in part to details of computation. In the earlier publication 
the curve of fertility m(a) was smoothed by the method of translation, with a Gaussian 
distribution as basis. In the method here presented the raw data were used without 
smoothing, except such as is inherent in the process of the calculation described. 

11 At any rate, Qio + Qi: must be small, since Q; + ... + Q, = 1.00313, and according to 
(24), with the convention that B(O) = 1, the sum of all the Q; must be equal to unity. 
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As a check, in order to be assured that no serious error was introduced in neg- 
lecting the remainder after dividing out the factors (y — y:) up to (y — y:), 
the product [Tia (y -— y) was computed and, after multiplying by a factor to 
make the absolute terms agree (.16635), was compared with the polynomial of 
(10). As a further indication, the coefficients of the product II were “decumu- 
lated” to obtain values of coefficients of the corresponding polynomial in z, to 


TABLE 4 


11 
Coefficients of Powers of y in Equation (10) and in the Product Ily - y3); 
1 
Also Coefficients of Powers of x in Equation (5) 


1 : " | Coefficients of x* in Equation (5) 
Coefficients of y | Found by Decumulation 


rs 
| 


In Equation (10)| In I(y— yi) | Of Column (2) | Of Column (3) 


| Se 
@) (3) | (4) 
| 
| 
| 


| 
| 


| 
| 
| 
| 


(5) 


— 1.00000 — .99915 
6.64127 6.64072 +.00014* | +.00065 
16.64550 16.64782 — .00057* | —.00432 
24 .34106 24.24197 | .02118* | .02398 


23 .16840 23.18070 .21781 21774 
15.05650 15.07338 | 33400 | .33354 
6.7250 | 6.73812 | 27427 | 27474 
1.99717 | 2.00316 | .18963 . 18882 
36404 36555 .10607 | 10641 

| 03483 | 03501 | 02268 .02276 
| 00127 00128 00116 00117 


| .00001 .00001 .00001 .00001 


. 16635 . 16635 


' | 


* In computing the denominator of Q according to (22) the values of the coef- 
ficients c, and c2 were arbitrarily made zero and the value of c; (age 15) was ad- 
11 


justed to .02073 to retain the total >) c; = 1.16635. 
1 


compare with values of c.. The results are shown in Table 4. In view of the 
fact that the (numerically) highest roots were determined only in first approxi- 
mation, the agreement is satisfactory. 

It is to be noted that instead of applying the solution (6) to compute values of 
B(t), these latter can, of course, also be obtained directly, by carrying forward 
step by step the original recurrent series; or, alternatively, the births in suc- 
cessive generations can be computed step by step and the total births obtained 
by addition. The advantage of the solution (6) is that it enables one, if desired, 
to obtain B(t) for any value of ¢ without having to compute B(é) for all inter- 
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vening values of ¢; also, the solution in an exponential series gives a better ideg 
of the general nature of the process, as well as a direct indication of its asympto- 
tic course for large values of ¢, when the first term Qoro ‘ with the positive real 
root 2» dominates all others. However this may be, it is interesting to compare 


TABLE 5 


Synopsis of Results of Computation of B(t) as 2 Qz-, Column (8), and as 
2B, (t), Column (9), where B,(t) = Births per Unit of Time in nth 
Generation at Time t. (Time Unit = 5 years) 


AnQrtaztl A oan Ss aac iain 
G G?+ H? } 


z* or u# .97322* 


17,716 
18,204 
18,704 
19 219 
19,748 
20,291 
20 ,850 
21 ,423 
22,013 
22,619 
23 ,241 
23 ,880 
24 538 
25 ,213 
25 ,907 
26 ,619 
27 ,352 
28 ,105 
28 ,878 
29 ,673 
30,489 
31,328 
32,191 





33 ,076 








— .19800 # |— .44720%| —.47587 # ZA 


28 ,482 
—415 
—19,498 
—15,222 
1,022 
11,057 
8,102 
—1,001 
—6, 248) 
—4,294 
792) 
3,519) 
2,265) 

— 568) 
—1,976 
—1,188) 
385 
1,106 
620) 
—251| 
—617 
—321 
160) 








343) 





24,234) 100,313 
—13,781 527 
3,329}  —274 
2,256} 2,326 
—3,362} 21,588 
2,223) 33,460 
—749| 27,470 
—169} 19,745 
445) 16,823 
| — 844) 18,012 
—194 | 145| 24,028 
—45 | —1) 27,328 
79 | —55| 26,841 
18 | 51| 24,706 
—32 | —26 23,878 
—8 | 25,424 
13 | | 27,757 
3 | 29 , 206 
—5 | | 29,498 
wa 
2 | 

1 
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TABLE 5—Continued 


B,(.) 





Generalions, n 


@ | ® | ® | © 





| 


ao) | ait) | @2 | 3) (14) 
cael | eeiinaiinbsleiita 





(9) | 


100,000 | 
0 | 
0 
2,072 | 
21,781 
33,398 | | 
27 ,472 | 43) 
19,866 £03) 
16,735 | 6,128 
| 17,954 | 15,685, 1 
24 ,033 23 ,889 28) 
27 ,361 | 27 ,022| 338 
26,878 | 24,905 1,973 | 
24 ,696 18,481) 6,214 1 
| 23,851 | 10,980 12,858, 13 
25,410 5,345 19,941, 124 
27,759 2,050, 25,030, 679 
29 ,219 526, 26,316 2,377 | | 
29 506 76 28,527) 5,897 | 6 
29,414 5| 18,092) 11,271 | 46 
29 ,862 12,041) 17,579 | 242 
31,000 6,906, 23,191 | 903 
32,348 3,381) 26,442 | 2,523 2 
33 ,423 1,397) 26,426 5,583 17 





OONIQarhwndDs KH © 




















the result of the computation by means of the exponential series, carried out as 
set forth above, with the corresponding results of the computation of births in 
successive generations. This comparison is exhibited in Table 5. 

It will be seen that the agreement is good except for the second to fourth 
items, where perhaps the omission of the terms contributed by the numerically 
highest roots makes itself felt. 


5. Discussion. 


(a) The real roots of the characteristic equation (5). It can be shown [8] that only 
one of the real roots for x can be positive, and that the absolute value of any 
other root must be greater than the positive real root. 
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The negative real roots which make their appearance in the numerical 
example call for special comment. Practically, the “higher” negative roots are 
of little importance, at any rate in this example—first because the constants Q 
with which they are associated are relatively small; second because large absolute 
values of negative roots imply rapid damping, so that corresponding terms Qr‘ 
very soon become negligible as ¢ increases. Thirdly, the determination of these 
roots would be subject to a wide range of uncertainty, corresponding to the large 
percentage fluctuations or errors of determination of the values of the functions 
p(a)m(a) = cq at the upper end of the reproductive period. 

But in theory these negative real roots suggest some pertinent questions. 
One wonders what would happen to them if the data were given, say, for single 
years of age, instead of 5-year groups. Instead of an equation of eleventh degree 
we would then have one of 55th degree. Furthermore, in those cases in which it 
may be permissible to pass to the limit, so that an integral equation takes the 
place of (2), negative roots for x would seem to be exciuded as they would make 
the integral in (13) meaningless. 

A problem of perhaps little practical importance but of some theoretical in- 
terest may arise here, to which reference has also been made by P.H. Leslie in a 
recent article in Biometrika,” in connection with a different procedure. 


(b) Effect of finer subdivisions of histogram of p(a)m(a). The effect of this on 
equation (5) for x is not obvious at sight, since new coefficients would be in- 
serted between previous terms. The effect is more easily understood from a con- 
sideration of equation (8) for y. Here finer subdivisions would introduce new 
terms only beyond the last term originally present. The original terms would 
not be changed at all in form, and those involving only lower moments would 
be changed but little in numerical value, provided that the original histogram were 
not so coarse as to give inappropriate values even for these lower moments. 
The result, then, of finer subdivision of the histogram, would be to change the 
computed values of the lower roots only in minor degree. But the four negative 
real roots, depending in considerable measure on the higher terms of (5) or (8), 
would presumably be materially altered, and might perhaps give place to further 
complex roots. In any case they would be followed by new roots even more 
remote from practical significance than the original eleven. 





(c) The result as an interpolation formula. Strictly speaking, the solution (6) 
of (2) is applicable only for integral values of ¢. In particular, terms arising 
out of the negative real roots of (5) for x are obviously not adapted to furnish 
interpolated values of B(t) for fractional values of ¢, since fractional powers of 

12 See [9] and [10]. For a brief summary and analysis of Leslie’s paper [9] see a review 
signed with the initials WGB in the Jourl Inst. of Actuaries Student’s Soc., Vol. 4 (1946), 
Part II. The first application of the matrix method to these problems seems to be due to 
H. Bernardelli, ‘‘Population Waves,’ Jour. of Burma Research Soc., Vol. 31 (1941), Part I, 
p. 1. 














erica] 
ts are 
nts Q 
solute 
| Q rt 
these 
large 
‘tions 


ions, 
ingle 
gree 
ch it 
: the 
nake 


l in- 
ina 


$s on 
. in- 
2On- 
new 
yuld 
uld 
yere 


the 
ive 
(8), 
her 
ore 











APPLICATION OF RECURRENT SERIES IN RENEWAL THEORY 205 





negative quantities in general are complex. Over the range of ¢ where the first 
real root together with the three parts of complex roots adequately describe the 
process under discussion, these terms alone are, in this sense and to this extent, 
suitable for interpolation, disregarding the terms corresponding to the other nega- 
tive roots. 

Even less suitable for interpolation purposes, it would seem, would be terms 
arising from further negative roots that might be introduced by a finer sub- 
division of the histogram of original data. If we suppose this subdivision carried 
to great lengths, and if negative roots still appeared under these circumstances, 
they would give rise to rapidly oscillating positive and negative terms for even 
and odd integral values of ¢ respectively (the time unit now being a subdivision 
of the original time unit) with no appropriate interpolation between these 
integral values. 

One further point calls for comment. In the process of idealization of the 
problem discussed, it has been assumed that p(a) and m(a) are independent of 
time, and the conclusions reached must be construed in the light of this assump- 
tion. In itself this would hardly call for comment, as it is a matter of common 
understanding. But the question does arise whether the assumption itself is 
free from implied internal contradictions. 

In a recent publication, P. K. Whelpton™ has drawn attention to the fact that 
in times of rapid changes in the birth rate, the assumption of age specific fertilities 
being held constant at the values observed in a given calendar year may imply 
that some of the women had more than one first child, a logical impossibility. 

The data used in the present numerical example are derived from a period of 
relatively undisturbed birth rate (1920), and do not involve any such conflict. 
But, in the light of Whelpton’s contribution one may ask the broader question 
whether the computation of an intrinsic rate of natural increase and related 
parameters based on age specific fertility as observed in one calendar year 
retain any practical value at all. 

In answering this question, two considerations will be weighed. First, that 
ordinarily the rates computed in the usual way differ but little from those ob- 
tained by taking into account order of birth as in Whelpton’s procedure. Sec- 
ondly, that the computation using over-all values of m(a) for all orders of birth 
combined is a relatively simple matter based on data commonly available; 
whereas the more complete treatment of the problem taking into account order 
of birth is considerably more complicated and often not possible at all for lack 
of detailed data. 
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SOLUTION OF EQUATIONS BY INTERPOLATION 


W. M. Kincaip 
University of Michigan 


Introduction and summary. The present paper deals with the numerical 
solution of equations by the combined use of Newton’s method and inverse in- 
terpolation. In Part I the case of one equation in one unknown is discussed. 
The methods described here were developed by Aitken [1] and Neville [2], but 
do not seem as widely known as they should be, perhaps because the original 
papers are not readily available. (A short summary of Aitken’s work will be 
found in a recent paper by Womersley [3].) Mention should also be made of 
an interesting paper by Spoerl [4], which treats the same problem from a some- 
what different viewpoint. 

In Part II these methods are extended to sets of simultaneous equations. 


PART I. EQUATIONS IN ONE UNKNOWN 


1. Nature of the problem. We first consider the problem of locating, to any 
desired degree of accuracy, a real root 2 of an equation of the form 


(1) y(x) = 0 


where y(x) is assumed to be analytic in an interval containing the root in ques- 
tion. Since we shall not be concerned here with the necessary preliminary work 
of separating the roots, etc., we may suppose that x is known to lie within a 
given interval that contains no zeros of y’(x). (Multiple roots are thus ex- 
cluded; but of course any such root is a simple root of an equation obtained from 
(1) by differentiation, and the methods described below can be applied to this 
equation.) 


2. Aitken’s method of interpolation. The method to be described, which 
may be regarded as a generalization of Newton’s, depends on the use of inverse 
interpolation. It is therefore desirable to recall a few points from the theory of 
interpolation before proceeding further. 

Let f be a function such that f(é) is known for ¢ = 4, &2,---,t,. Then the 
Lagrange interpolating polynomial fi»...,(/) is defined by 





. +-96< 8 *(t — th) 
fia) = §0) (yg = &) (= be) 


+ fb (i — “a: — = . — tn) 


as 
ot 
207 
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We note that 
f(t) t — &| 


i. t—b nvt—-t 'f(bt—b 

a(t Se. secs aN reared 
fiall) = fh) ——F + fle) =F ae 
fe(h ty | fie3---n—1(t) t = ty | 
fo,(t) t — t . . 23---n(b t— é, | 
— (i) t - : - ’ <= ’ fir..-n(t) = | fi — ( ) a ee, 

ty = ts ty _, tn 
so that fis3....(f) can be evaluated for any given value ft of ¢ by a succession of 
linear interpolations. It is convenient to arrange the work in a table like the 
following (n = 4): 


fi 23 (t) 


TABLE la 
f(t) I II | Parts 


— - —$—$—$— $$ 


f(t) to-t 
Fia(to) 
F(t) Fiza(to) tote 
fo3(to) 
ls F (ts) Fosa(to) tots 
tsa(to) 
ty S(t) lols 
This form is well adapted for machine computation, for each denominator 
t; — t; = (t — t;) — (t — ¢t;) automatically appears in one set of counters when 
the corresponding numerator is obtained in the other. 
If f’(t) is known at one or more of the given points, this information is readily 
fitted into the scheme. For we see that 


(4) fu) = lim f(t) = f(a) + @ — Of’) 


tot, 
and all that is necessary is to repeat certain entries in Table Ia and to fill in col- 
umn I by using (4) as indicated in Tabie Ib. The extension to higher deriva- 
tives is obvious. 
TABLE Ib 


II | Parts 


to-t, 
Fir(to) 
firo(to) | to-h 
fio(to) | fires(to) 
fi2s(to) to-te 
fos(to) fi233(to) 
foss(to) lo-ls 
fs3(to) 
to-ts 
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In applying the above to obtaining the root x of (1), we must suppose that 
y(z) is tabulated or can be computed for a set of values of x in the neighborhood 
of %. What we do not know is the value of x corresponding to y = 0. It is 
therefore convenient to regard x as a function of y whose value is known at cer- 
tain points and then interpolate to get x = x(0). That is, we let y take the place 
of tand x that of f(t) in the preceding discussion, while 0 replaces 4. The work 
is slightly simplified by the fact that the column of “parts” becomes identical 
with the left-hand column which contains the y’s and can therefore be omitted. 


3. Application to an example. The procedure will be most clearly indicated 
by an example. Consider the equation 


(5) y = 2+ 22° — 52? — 8x +1 = 0, 


which has a root between 0 and 1. (If the root were located elsewhere, it would 
be desirable to shift it to this interval in order to simplify the computation of y.) 

The work of evaluating this root to ten places is summarized in Table II, and 
explained below. In the first column, the numbers in parentheses are values 


of oe and the other numbers are values of y, corresponding to the values of x 


in the second column. 
TABLE II 


Vv x I 


1.000 000 000 000 0.000 | 
(—8.000 000 000) 0.000 | 0.125 009 000 00 | 
0.152 100 000 000 0.100 | 0.117 938 436 13 | 0.116 671 702 00 
—0.001 054 385 279 | 0.117 6 882 964 17 | 884 075 87 
(—9.081 459 548) 0.117 3 896 94 | 3 890 52 | 0.116 883 877 O1 
0.008 022 855 936 | 0.116 3 842 98 | 90 67 90 68 
(—9.073 020 416) | 0.116 4 254 15 | 90 74 90 68 


Xo = 0.116 883 890 7 


The procedure is as follows. Taking x = 0 as a first approximation to 1, 
we find that y(0) = 1, y’(0) = —8, and record this data in the y and x columns 
of the table. Note that for convenience, the value of y’(0) takes the place of a 
blank entry in Table Ib. We now apply (4), which here takes the form 
(6) an(0) = 2 v= + (0 -— 1) i La =0+ a =0+ = = 0.125 

dz z=0 
and enter the result in column I. Note that this is equivalent to one step of 
Newton’s method. 

In view of (6) we take x = 0.1 for our next approximation and apply (3) to 
obtain the second entry in column I and the first in column II. This last sug- 
gests « = 0.117 for our next trial value. (We do not compute y’(0.1), as little 
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would be gained by doing so, and the time is better spent in going ahead ag 
indicated.) Finding y(0.117) and filling in the table gives us the root to six 
places. 





4. Employment of tables. Continuing in the same line, it would seem nat- 
ural to take x = 0.116884 at the next step; and doing so would lead to the most 
rapid convergence. But another consideration enters. Up to this point the 
values of y were computed with the aid of the WPA Table of Powers, which igs 
limited to three places in the argument. Rather than going to the extra labor 
of evaluating y(.116884), we proceed as indicated in the table, using y’(.117), 
y(.116) and y’(.116), and stopping when the values of zx in the last column 
agree to the desired number of places. 

This point has been dwelt on because it is likely to arise whenever tables are 
used in evaluating y(x). In the example just given, to be sure, we had a certain 
freedom of choice; but if y(x) is not algebraic, direct computation may be quite 
impractical. It may be noted that in such cases the method of inverse inter- 
polation is not only faster than the simple Newton’s method but is capable of 
giving more accurate results. 

The error in the final result can be estimated from the standard formula for 
the error of interpolation, but this may be awkward because it requires the 
evaluation of higher derivatives of x with respect toy. In practice it is generally 
safe to rely on agreement of different interpolated values, and of course the result 
may be checked by substitution in the original equation. One simple point is 
worth noting, however—if the error in the original column of z’s is O(e), that in 
the successive columns to the right is O(é), O(e), ete. 





















5. Applicability of the method. Although the example we have presented is 
algebraic, the method is, of course, equally applicable to transcendental equa- 
tions. Moreover, it can be used, theoretically at least, to yield complex as well 
as real roots. The sole difficulty is that the numerical work becomes cumber- 
some in this case; how serious it is depends on the type of computing machines 
used. If the equation is algebraic, Bernoulli’s [5], [6] and Graeffe’s [7] methods 
are applicable. In fact, they are likely to be the most effective since they do not 
require prior knowledge of a first approximation to the root. If the alternative 
procedure of replacing the equation by two simultaneous equations for the real 
and imaginary parts of the root is decided upon, the methods described in the 
next section may prove useful. 


PART II. SETS OF SIMULTANEOUS EQUATIONS 


6. Two equations; general considerations. It is natural to take up next the 
problem of finding the simultaneous solutions of two equations in two unknowns. 
Let these equations be 


(7) 7 u(x, y) = 0, v(z, y) = 0, 
where u and v are analytic functions of x and y. 
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If we had a general method of interpolation of functions of two independent 
yariables, the problem could be solved in a fashion similar to that used in the 
preceding section. That is, u and v would be computed for values of x and y 
near the desired ones; then x and y would be regarded as functions of u and v 
and interpolations would be performed to obtain the values corresponding to 
@=v = 0. 

It is easy to set up interpolating functions in a variety of ways, but the author 
has found none that are satisfactory for the problem in hand. Note that what is 
required is to determine the value of a function at any point in the plane, given 
its values at a set of fixed points. The most obvious idea is to use polynomials of 
the least possible degree for this purpose, as is done in the case of a single variable. 
In this case, however, the coefficients of a polynomial of the nth degree are de- 


termined by its values not at n + 1 but at i points; thus if a func- 


tion is given at 5 points, no unique quadratic interpolating polynomial can be 
constructed. What is worse, even if a function is given at 6 points, say, the 
quadratic polynomial determined will in general have large coefficients and take 
on unreasonable values if all the points happen to lie close to a common conic. 
Other schemes considered by the author have similar drawbacks, though the 
possibility of course remains of finding a suitable one by further research. 

The problem can also be handled; at least in principle, by eliminating one of the 
variables; but, apart from the difficulty of carrying this out in practice, the 
resulting single equation is likely to be more complicated in form than the original 
two. If so, solving it may require more computation than would be involved in 
attacking the original equations directly by the methods described below. So 
far is this true that even when a single equation is given in the first place it may 
be advantageous to replace it by a set of simpler equations. 


7. Newton’s Method. Although a direct extension of the method of inverse 
interpolation is not presently available, Newton’s method may be suitably 
generalized for this case. 

Starting with equations (7), we set up the auxiliary variables 


(8) X = Uy — VW, Y = w, — vuz, 


where the subscripts denote partial derivatives; uz = = » Wy = a , etc. 


We have 


a — vu _ 
dz 2vy zUy zy zy» ay vy UUyy ; 
(9) 

oY oY 


—_ = Uvez — VUzzr ’ —— 
oy 


5 = Uy Ve — Vy Us + Way — Uy. 
x 
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For u = 0,v = 0, equations (8), (9) reduce to 


oy Ox Ox oy , 


where J is the Jacobian of u and v with respect to z and y. 

Equations (10) will hold approximately for values of x and y near those satis. 
fying equations (7). That is, in the neighborhood of a solution X can be re- 
garded as a function of x alone and Y as a function of y alone. Then if z = x, 
y = Yo is the desired solution, (1 , y:) is a point in its neighborhood, and z, = 
X (x ? Yi); Y, oe Y (x1 ? Yi); Ji = J (x ? Yr); we have 


4 Y 
(11) t9™~%1 = - ) Yo™ Yi +—. 
Ji Ji 


Also if (a2 , y2) is another point near (x, yo), 


(12) i at X2 — %2 Xi 4": Y2 — yo Yr 


a ip Oe ene 
X 2 = X1 7 oe Y; 


Relations (11) and (12) can be used to obtain successive approximations to the 
solution. Use of these relations corresponds to employing Newton’s method 
and linear interpolation for the solution of one equation in one unknown. 

As a first example we consider the equations 


usxtaty-— 
13) 
( (edt ois 


We have 


Ue=2re+y Uy = x + 2y 
(14) ; 
vz, = Qry Vy = x + 2y. 


Drawing a rough graph indicates a solution near (1, 1). We evaluate u, », etc., 
at this point as shown in Table III. Using (11) we get (2, 0) for a second approxi- 
mation, and proceed as before. We can now use both (12) and (11) to get new 
approximations; they are (1.33, 0.57) and (1.25, 0.50), and are entered in the 
last two columns of the table. We therefore try (1.3, 0.5) next, and continue 
in this fashion until the desired accuracy is attained. Both (11) and (12) are 
used at each step and the values obtained are entered in the last two columns. 
The entries in the numbered rows are obtained by using (11), the others by using 
(12). The number of places to take in each succeeding step is judged from the 
agreement shown. 

Table IV indicates the process of finding a second solution of (11) by the same 
method. The convergence is very rapid in this case, mainly because the first 
guess is fairly close. 
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8. Inverse interpolation. In the preceding section, attention was drawn to 
the difficulty that may arise when tables, necessarily limited to a certain number 
of places in the argument, are used in the computation. In the example just 
discussed the values of u and v were easily computed directly to the number of 
places wanted. But a glance at the work will show that if we had been limited in 
computing wu and v to values of x and y having, say, two decimal places, the solu- 
tions could have been carried to four places only. 

The device adopted in the preceding section was to use quadratic and cubic 
interpolates to secure greater accuracy, and it might occur to us to try the same 
idea here. But for such an interpolation to be strictly valid, equations (10) 
would have to hold identically. Since they hold only approximately, an error 
is introduced which, in general, is of the same order of magnitude as the error in 
linear interpolation. Thus continuing the interpolation would not improve the 
results. 

However, this very situation suggests a way out. For suppose we give x a 
constant value x; , and compute X and Y for a number of values of y. Forz = 
a,, both X and Y can be regarded as functions of y alone; or we can regard X 
and y as functions of Y. Doing so, we can interpolate to any number of stages 
to find values of X and y corresponding to Y = 0; call these X,,y,. Assigning 
z other constant values x2, 13, °°: , Um, Wwe repeat the process, getting a set 
of values X2,--- Xm and ye,--*, Ym, all corresponding to Y = 0. Now 
along the curve Y = O we can regard x and y as functions of X; performing one 
more interpolation, we obtain the desired values of x, y corresponding to X = 
Y = 0. The error in the final result can be estimated from the errors in the 
interpolations, and is of the same order of magnitude as the greatest of these. 

It will be noted that we did not refer to the definitions of X and Y in describing 
this procedure. Any pair of independent (analytic) functions X’ and Y’ having 
the property that X’ = Y’ = 0 when u = v = O could be used. However, it 

/ / 
is convenient to choose them so that - and _ are small. Probably the 


simplest course is to set 
X’ = au + by, Y’ = agu + bo, 
where a; , a , b; , be are constants such that 


ay, vy ae Uz 
~S = —_— 


by Uy be Uz 


Let us apply this procedure to the example we have already worked (Table III). 
Suppose we wish to use values of x and y having not more than two decimal 
places. Within this restriction, we can still carry through the first few steps 
indicated in Table III to ascertain that 1 ~ 1.514, yo ~ 0.375 where (2 , yo) 
is the desired solution. At the point (1.51, 0.37) we have 


X = 3.020lu — 2.25v, Y = 1.1174u — 3.39v. 
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Noting the ratios of the coefficients of «u and v, we select 
X’ = 4u — 3v, Y’ = u — Bv. 


Next we evaluate X’ and Y’ for the 16 points having x-coordinates 1,50 
1.51, 1.52, 1.53 and y-coordinates 0.36, 0.37, 0.38, 0.39, as shown in Table v. 
Starting with the four points for which x = 1.50, we interpolate to find the values 
of y and X’ corresponding to Y’ = 0; they are y,; = .3750000007, X; = —.1406- 
250025. We proceed in the same way with the points corresponding to the other 
values of x; the results, as shown, are y2 = .3749981706, - C= — 03908741039. 
y3 = .3749927660, X3 = .06302572977; ys = .3749839124, X, = .1657149545, 
(The extra digits given in Table V are to take care of rounding-off.) Finally, 
using these values, we interpolate to find the values of x« and y corresponding to 
X’ = 0, and get 


x = 1.5138345192, y = .3749965140. 


Comparing these results with those obtained earlier, we see that they are in error 
by about 1 unit in the ninth place; a distinct improvement over the four correct 
places that could have been secured without using this device. Note that if we 
had not had our earlier results for comparison, a check could have been obtained 
by carrying through the interpolation in the reverse order; i.e., starting with 
fixed values of y and finding values of x and Y’ corresponding to X’ = 0. 

As in the case of one equation in one unknown, derivatives could be brought 
into the interpolation scheme, permitting greater accuracy with fewer points. 

Ox O02 07x 


But the derivatives needed would be aX’? ay’ aX’ay’ etc., and the general 


setup would be rather awkward, so that extra labor would probably be required. 


9. Three or more equations. The methods discussed in this section are 
readily extended to the solution of three or more simultaneous equations in an 
equal number of unknowns. For example, if we are given three equations of the 
form 


u(x, y,2z) = 0, v(z,y,z) = 0, wa, y,z) = 0, 
we define new variables 
uvw Uz Uz Wz 
X =| Uy% WW, |, : luvw P Z | Uy Vy Uy 
Uz Vz Wz Uz, Vz Uz uv WwW 


which are analogous to the X and Y of (8); from this point on the work is practi- 
cally the same as before. 


REFERENCES 


{1] A. C. Arrken, “On interpolation by iteration of proportional parts, without the use 
of differences,’’ Edinb. Math. Soc. Proc., ser. 2, Vol. 3 (1932), pp. 56-76. 





SOLUTION OF EQUATIONS 219 


(2) E, H. Nevitve, “Iterative interpolation,”’ Indian Math. Soc. Jour., Vol. 20 (1934), pp. 
87-120. 

(3] J. R. WoMERSLEy, “Scientific computing in Great Britain,’ Mathematical Tables and 
aids to Computation, Vol. 2 (1946), pp. 110-117. 

[4] C. A. SPOERL, ‘‘Solving equations in the machine age,’’ Amer. Inst. Actuar. Record, 
Vol. 31, Part I (1942), pp. 129-149. . 

(5) T. C. Fry, ‘Some numerical methods for locating roots of polynomials,” Quart. Appl. 
Math., Vol. 3 (1945), pp. 89-105. 

(6] W. M. Kincarp, ‘‘Numerical methods for finding roots and vectors of matrices,’’ Quart. 
Appl. Math., Vol. 5 (1947), pp. 320-345. 

\7] E. Bopewie, ‘On Graeffe’s method for solving algebraic equations,’’ Quart. Appl. 
Math. Vol. 4 (1946), pp. 177-190. 





/ 
4 


ESTIMATION OF A PARAMETER WHEN THE NUMBER OF 
UNKNOWN PARAMETERS INCREASES INDEFINITELY WITH 
THE NUMBER OF OBSERVATIONS 


By ABRAHAM WALD 


Columbia University 


Summary. Necessary and sufficient conditions are given for the existence. 
of a uniformly consistent estimate of an unknown parameter @ when the succes- 
sive observations are not necessarily independent and the number of unknown 
parameters involved in the joint distribution of the observations increases in- 
definitely with the number of observations. In analogy with R. A. Fisher’s 
information function, the amount of information contained in the first n observa- 
tions regarding @ is defined. A sufficient condition for the non-existence of a 
uniformly consistent estimate of @ is given in section 3 in terms of the information 
function. Section 4 gives a simplified expression for the amount of information 
when the successive observations are independent. 

2. Introduction. J. Neyman has recently treated the following estimation 
problem’: Let X,, X2,--- , etc. be a sequence of independent chance variables 
the distribution of each of which depends on some unknown parameters. Two 
kinds of parameters are distinguished, structural and incidental parameters. A 
parameter @ is called structural if there exists an infinite subsequence of the 
sequence {X,} such that the distribution of each of the chance variables in the 
subsequence depends on @. Any parameter which is not structural is called 
incidental. Neyman has considered the case when there are a finite number of 
structural parameters, say 6,,--- , 0, and an infinite sequence {£:}, (¢ = 1, 2, 

- ,ad inf.), of incidental parameters. He has studied the problem of consistent 
and efficient estimation of the structural parameters and has obtained several 
interesting results. He has shown, among others, that the maximum likelihood 
estimate of a structural parameter @ need not be consistent, even when consistent 
estimates of @ exist. Neyman has also given a method for obtaining consistent 
estimates of the structural parameters. This method, however, is applicable 
only under certain restrictive conditions. 

In this paper we shall consider a more general case than that treated by Ney- 
man, but we shall concentrate on one aspect of the problem, namely that of the 
existence of consistent estimates. 

Let {X.}, (@ = 1, 2,---, ad inf.), be a sequence of chance variables, not 
necessarily independent of each other. It is assumed that for each n the chance 
variables X,,---, X, admit a joint probability density function 
Pn(t1,°°* ,%n| 0,8, +--+ , &n) where 0, & , &,--- , etc. are unknown parameters.’ 


1 Address given by J. Neyman at the meeting of the Institute of Mathematical Statistics 
in Atlantic City, January, 1947. 
2 While @ is assumed to be a real variable, we admit é; to be a finite dimensional vector, 
i.e., & = (1, ..., E:&;) where k; may be any finite positive integer. 
220 
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We shall require that the consistency relations among the density functions 
fh, p2,°** ,etc. be fulfilled, i-e., 


+00 
(1.1) / Pati Ain = Da, (n = 1,2, ---, ad inf.). 


It should be remarked that it is not postulated that p, actually depends on all the 
parameters that appear as arguments in p,. It is merely assumed that pp, 
does not depend on any parameter that does not appear as an argument in p, , 
i.e., Pn does not depend on é; for anyz > n. It follows, however, from (1.1) that 
if p, depends on a parameter &, then also p,», depends on é for any m > n. 

Neyman’s definition of structural and incidental parameters can be extended 
to the case of dependent observations considered here by saying that the dis- 
tribution of X; does not depend on a parameter é if and only if the conditional 
distribution of X; for any given values of X,, --- , X;_; does not depend on &. 
It is not postulated that each of the parameters £, , &,--- , etc. is incidental; 
some of them may be structural. We shall not make an explicit distinction 
between structural and incidental parameters, since for the purposes of the 
present paper this does not seem to be necessary. 

In this paper we shall deal with the problem of formulating conditions under 
which a uniformly consistent estimate of 6 exists. A statistic t,(a@1,--- , 2n) is 
said to be a uniformly consistent estimate of 6 if for any positive 6 
(1.2) lim prob. {|t, — 6] < 6} =1 
uniformly in @ and the é’s. 

In section 2 a necessary and sufficient condition is given for the existence of a 
uniformly consistent estimate of 6. In section 3 the amount of information 
supplied by the first n observations concerning @ is defined. It is then shown 
that if the amount of information is a bounded function of n over a non-degener- 
ate 0-interval, no uniformly consistent estimate of @ exists. Section 4 gives a 
simplified formula for the amount of information in the case when the X’s are 
independently distributed. 

2. A necessary and sufficient condition for the existence of a uniformly 
consistent estimate of 6. In deriving a necessary and sufficient condition for 
the existence of a uniformly consistent estimate of 6, use will be made of some 
results contained in a publication of the author [1] dealing with statistical decision 
functions which minimize the maximum risk. In [1] it is assumed that the 
domain of each of the unknown parameters is a closed and bounded set and that 
Pn is continuous jointly in all of its arguments. Thus, in order to be able to use 
the results obtained in [1], we shall have to make the same assumptions here. 
In what follows we shall, therefore, assume that each of the parameters @, & , 
f,---+ , etc. is restricted to a finite closed interval and that p, is a continuous 
function of 271, °-- , 2n, 9, £1,°°*, &n. 

Let [a, b] (a < b) be the 6-interval to which the values of @ are restricted. 
Clearly, if tn(a1,--°- , 2n), (n = 1, 2,---, ad inf.), is a uniformly consistent 
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estimate of @, then also t* isa uniformly consistent estimate of @ when i. = th 
when a < t, < b, (* = awhent, < aandt, = b whent, > b. Thus, without 
loss of generality, we can restrict ourselves to estimates t, which can take values 
only in the interval [a, b]. Uniform consistency of ¢, is then equivalent with the 
condition 

(2.1) lim El(t, — 0)°|0,&,---,&] = 0 

uniformly in @ and the é’s. For any chance variable u the symbol 
E(u | 0, &: ,€,--°-) denotes the expected value of uw when @, & , &,--- are the 
true parameter values. 

In [1] a non-negative function W(t, , 0), called weight function, is introduced 
which expresses the loss suffered when ¢, is the value of the estimate and @ is the 
true value of the parameter. The risk is defined in [1] as the expected value of 
the loss, i.e., the risk is given by 


(2.2) rn(8, & i. En) = E(W (t, ’ 6) | 6, £1 Sores &,]. 
If we put W(t, , 0) = (tn — 0)”, we have 
(2.3) rn(0, & aie En) =— El (tn os @)” | 0, & en. ae &,]. 


It can easily be verified that Assumptions 1-4 in section 3 of [1] are fulfilled 
for the weight function W(t, , 0) = (tn — 0). Thus, all results obtained in 
[1] can be applied to the risk function given in (2.3). According to Theorem 4.1 
in [1] the risk function given in (2.3) is a continuous function of 0, &,--- , & 
for any arbitrary estimate ¢,. We shall denote the maximum of (2.3) with re- 
spect to 6, &,---, & by r,[t,]. Thus 1r,[t,] is a functional which associates a 
non-negative value with any estimate function ¢, . 

It follows from (2.1) that ¢, is a uniformly consistent estimate of 6 if and only 
if 
(2.4) lim r,[t,| = 0. 

For any @ and for any n let F,.(&,--- , 2 | @) be a cumulative distribution 
function of &,--:,&,. Let, furthermore, 


Gn(X1 ee | 0, F,) 


(2.5) + +o 
=| cee [ Dalti, +++, Xn} 0, bi +++ En) P(E, -- , 


— 6 ia 


We do not require that F,, F2,--- , etc. satisfy the consistency relations, i.e., 
lim Frii(fi, -**, €n41| 0) is not necessarily equal to F,(&,---, & | 6). 


En+1=0 

3 In verifying Assumption 4, we may assume that p, is always > 0, since for any given 
values 6, &},..., &: we may restrict the domain of (z,..., xn) to the subset of the sample 
space where p» > 0. 
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Hence, also the distributions g, do not necessarily satisfy the consistency rela- 
tions. Clearly 


° +0 +00 
(2.6) raltn] = [ sa [ (t, — 8)°qn(ars p°** ae | 6, F,,) day »°°* dry 


for any @ and any F,,. Hence, (2.4) and (2.6) imply that if ¢, is a uniformly 
consistent estimate of 6, then ¢, remains a uniformly consistent estimate of 6 
also when q, is the distribution of X,, --- , X, for any arbitrary choice of F, . 

For each n let C,(0, &, --- , 2) be a joint cumulative distribution function of 
0, &,°*:,&,. If this is regarded as an a priori distribution of 6, &,,--- , &, 
and if our aim is to choose ¢, so that 


E(t, — 6)” 


(2.7) - a — 
= [ | (t, — 0) prlti, --- ,%n|0,8,°++,&) dC, da, --- dz, 
is a minimum, then the best choice of ¢, is to put it equal to the a posteriori mean 
value of 6. Let tn(a,---,2%n3 Cn) denote the a posteriori mean value of 6 
when C,, is the a priori distribution, i.e., 


| ep ates, pyere ,tn| 0,8, aie ks , €n) dC, 


(2.8) th (ay 9. Ss ee C,,) — _ 

| Pal yt Un | 6, & ae én) dc, 
where the integration is to be taken over the whole domain of the parameters 
0, f&,-°-:,é&,. Let, furthermore, 7,[C,] denote the value of (2.7) when t, = 


tn(%i, °°: ,%n3;C,). According to Theorem 4.4 in [1] there exists a particular 
distribution C®, , called a least favorable distribution, such that 


(2.9) FalCa] S FAC] 
for all C,. Let 
(2.10) (t1,-°+,%n) = trl, --- 2,03 C%). 
It follows from Theorems (4.5) and (5.1) in [1] that for any estimate t, we have 
(2.11) ralt) 2 alle] = #.(C%). 


Hence, a necessary and sufficient condition for the existence of a uniformly 
consistent estimate of @ is that 


(2.12) lim 7,[C°] = 0. 

Let F.i(é&,---, & | 8) denote the conditional cumulative distribution of 
f,---, € for given @ that results from the joint distribution C,(0, & , --- &,) 
and let F%.(é ,-:+ , &. | @) correspond to C°(0, , ---, &). Clearly, any uni- 
formly consistent estimate of @ with respect to pr(a,--+ , tn| 0, &,°°: , En) 
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is a uniformly consistent estimate also with respect to gn(%:,-°-- , tn| 96, F,) 
for any F,. On the other hand, if qn(1,--- , tn | 0, Fx) admits a uniformly 
consistent estimate of 6, equation (2.12) must hold and, therefore, p,(2 , --- , 
XYn| 0, &,-°++, &) admits a uniformly consistent estimate of 6. Hence we 
arrive at the following theorem: 

THEOREM 2.1. A necessary and sufficient condition that 


Dn(t1, °°, tn | 0,8, °° 2) 


admit a uniformly consistent estimate of @ is that qn(ti1, --- , tn | 0, Fn) admit a 
uniformly consistent estimate of @ for any arbitrary choice of F, . 

3. Amount of information contained in the first n observations concerning 
the parameter §. We shall make the following assumptions: 

Assumption 1. The first two derivatives of p,(a1,-°--, tn| 0, £1, °°-, &) 
with respect to 6 exist. 

Assumption 2. We have 


+00 +00 ep 
(3.1) [ oes | Max | -5"| dx, --- dx, < @ 
Laas Lee 6 00 | 


and 


(3.2) | see Max ie _ (dx, +++ dt, < 0 
S Lo | OF 
for any n. 
Assumption 3. The integral 


- © log qn(ai, +++ 5 tn| 6, Fr 
[ nea | a ar 228918 Fd) dn(ti, «++ 2x} 0, Fa) dx, +++ di, 


exists for any 6, F,, and n where q, is defined by (2.5). 
Since 


ae dn on 


00 


A loggn 1°0°¢n c log =] 


and since, because of Assumptions 1 and 2, 


2 : ” On 
/ [SB diz, --* dx, = 0, 


we have 


© ° 3 lo In 
Lie meme 


+00 2 
9 log dn 
vs [ ( =a Gn At, +++ dry. 
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Let 


a f as ie log Gn 
(3.4) c,(@) = — \7 [. vee [. (2 ees) Qn dm --° az} 


Clearly c,.(@) = 0. We shall now show that 
(3.5) Cn41(9) = Cn(8) forn = 1,2, --- , ad inf. 


In fact, we can write 


= log qnsa(ti, +++» tn 18, Pass) _ _ OF log gn(mi, -*+ , tn 10, Fr) 
| og? 0g? 
3.6 
_ & log fasaCtnsi | tr, +++, En» 9, Faas) 
06? 





where Fr = lim Fay(éi,-++ , &n41 | @) and fn4r(@n41 | 1, °°*, Tn, 0, Fass) 


Ent+1=0 
is the conditional probability density function of X,4: given the values of 2, 
+++ , 2, and assuming that the joint density function of Xi, --- , Xa4 is given 
by Qn4i(t1,°°* »Xn4i| 0, Pn4i). Sincec,(6) S expected value of 


- a log qn(t1, +--+ , tn| 9, Fr) 
OF” 
2 
: -_ "t} is = O, inequality (3.5) must hold. 
In analogy with R. A. Fisher’s information function, we shall call c, (0) the 
amount of information contained in the first n observations regarding 0. We 
shall now prove the following theorem: 
THEOREM 3.1. Jf lim c,(6) S ¢ < © over a finite non-degenerate 0-interval I, 


and since the expected value of — 


then there is no uniformly consistent estimate of 0. 
Proor. If for any n, c,(@) S ¢c < © over the interval J, for each n there exists 
a distribution F,(&,--- , &. | @) such that 


ae on | als [ 8° log qu(ti, +++» tn 19, Fa) 
(3.7) fe 0 or 


- Gn(ti, +++, Xn|0, Fn) dm +--+: dz, Se+1 


for all n and for all @in J. Let ¢, be any estimate and let 


b,(8) = E(t, -)= | hi | (t, — O)qn(r1, coe Be | 6, F,) dx, --- dz, 


(3.8) “ ‘“ 
7 l oe l tngn(t1 , ++» tn 0, Pq) da +++ da, — 0. 


db,(6) 


Since ¢, is bounded, it follows from Assumptions 1 and 2 that 7 


exists and is 
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a continuous functionf{of 6. According to a theorem by Cramér [2] we have 


— ‘ (1 + Se) 
(3.9) E(t, — 0) = | cee | (t — 0)°an dx, -+- dx, = as _99 


= c+ 1 


for all 6 in J. Thus,gin order that lim E(t, — 6)” = 0 uniformly in 6, we must 


N= 


have 


. db, (8) 
3.10 ie 
“a — a 
uniformly in @ over J. Let J be the interval ranging from g toh (g < h). From 
(3.10) it follows that 


(3.11) lim [b,.(h) — bn(g)l = g — hk. 


TkaemOO 
Hence 


lim inf max [b,,(@)|’ = (g — hy 
n= @inIJ 4 
Since E(t, — 0) = [b,(0)]’, E(tn — 9)” cannot converge to zero uniformly in 6 
and Theorem 3.1 is proved. 

4. Formula for c,(6) when pn(ti,---, tn| 6, &,-°--, &) is equal to a(x 
| 0, £1) g2(X2 | 0, £2) ag Pn(Ln | 0, é,). Let gi (Xi Zip *** » Sea, & F,) be the 
conditional probability density of X; given 21, --- , «i. when the joint density 
function of 7, --- , Xn is given by gn(a1,---,an|0, Fn), (i Sn). Clearly, 


(4.1) -2(? log - od : E ee log ae) 


oe? 
Now 


(4.2) g(a; | M,°°* , 41,9, F,) = | gila; | 0, &:) dH(é; ime n °** 5 Mee, F, F,) 
where H;(é; | 21, °° , ti, 6, F,) denotes the conditional cumulative distribu- 
tion of &; given 2, --* , ti-1, assuming that F,.(& , --- ,&, | 0) is the joint cumu- 
lative distribution of & ,--- , &, and pr(ti,--: , @n| 6,61, --+ , &,) is the joint 
density of X,,--- , X, for any given values of 0, £1, --- ,&n,. 

It follows from (4.2) that 


+0 
a | 
-[ ee 9 9; da; = Cns(6) 


so] 0 "tog | on &) dC (&)  .« 
= g.1.b. | — [ ——— 9 _ . 9; dC; 


Cy (Ex) o— 0 
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where C;(&;) may be any cumulative distribution of &;. Hence 


lll a log «| - 
(4.3) gib.| B( age) | = Cas) 


and, therefore, 


(4.4) Cn(0) i a Cas(6). 


The quantity c,;(@) is simply the amount of information contained in the ith 
observation alone. Thus, formula (4.4) says that if X,,---, X, are inde- 
pendent, the total information contained in the first n observations is equal to 
the sum of the amounts of information contained in each of these observations 
singly. 
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INVERSION FORMULAE FOR THE DISTRIBUTION OF RATIOS 


By JoHN GURLAND 


University of California, Berkeley 


1. Summary. ‘The use of the repeated Cauchy principal value affords greater 
facility in the application of inversion formulae involving characteristic func. 
tions. Formula (2) below is especially useful in obtaining the inversion formula 
(1) for the distribution of the ratio of linear combinations of random variables 
which may be correlated. Formulae (1), (10), (12) generalize the special cases 
considered by Cramer [2], Curtiss [4], Geary [6], and are free of some restrictions 
they impose. The results are further generalized in section 6, where inversion 
formulae are given for the joint distribution of several ratios. In section 7, the 
joint distribution of several ratios of quadratic forms in random variables 
X,, X2,-°-- , X» having a multivariate normal distribution is considered. 


2. Introduction. We shall write 


f¢ ee f ola, t, -++t,) di dh --+ dt, 


= lim []-- | oe,e, +++, ty) dt; dt, «++ dt,, 
€;-0 


Timo eg<cl|[ti|<Ti 


which might be called the repeated Cauchy principal value of 


[| | gti, to, «+, tr) dt --- dt,, 


and which we shall use frequently. The results of this article may be regarded 
as extensions of the following theorem proved in section 4. 

THEoREM 1. Let X,, X2,°-:, Xn have the joint distribution function 
F(x, %2,°** , Ln) with corresponding characteristic function o(h, t,°-° , th). 
Let G(x) be the distribution function of (a,X1 + +--+ + anXn)/(b1X1 + --- + OnX), 
where a; , @2,°** ,@n,61,b2,-°-: , bn are real numbers. If 


P{d b;x; < o} = 0, 
1 
then 


(Q) Gl) + Ge -0) =1-= f Plier mi. a Han = bn2)h oy 


3. An inversion formula for distribution functions. Let /'(x) be a distribu- 
tion function and ¢(t) be the corresponding characteristic function. Then the 
following inversion formula holds: 

228 
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() F(t) + Fé — 0) =1- if eH gy 


PROOF. 


(f+ [)esoot= (fo +f) a [em are 
= [ dF (x) ([. + [) fr? :, 


by the Fubini theorem on the inversion of integrals. But 


& f pit(e-® . = sgn (x — &), 


md 


T . 
where sgn y = —1, 0,1 according asy <0,y=0,y>0. Since [ me dtis 


r t 
uniformly bounded in 7’, the principle of bounded convergence for Lebesgue 


integrals implies that 


i lim P dF (x) ([. + | ane : = [ sgn (x — &) dF(zx) 


Tt 440 Yo . 
(f- ° I. + [) sgn (x — £) dF(2) 


—F(é — 0) + 1 — F@). 
The required result follows at once. 


Another form of (2) may be obtained as follows: Let H(x), K(x) be distribu- 
tion functions, and y(t), x(t) the corresponding characteristic functions. Setting 


F = H,¢@ = vy, — = 0 in (2) yields 
Ho) + HO -0) = 1-14 wo, 
Tt t 
while setting F = K,¢ = x, in (2) yields 
K() + K(f —-0 =1- f x(t)e dt 
11 t 


Clearly 


_ =e 
@ K®+KE-0 = HO) +H0-0 +1 f OW xe PM 


If H = K, then y = x, and (3) reduces to a well-known inversion formula (cf. 
Kendall [7, p. 91}). 


4. Distribution of the ratio (a,X; + --- + anXn)/(byXi + +--+: + On Xn) 
with denominator positive. THrorem 1. Let X;, X2,---, Xn have the joint 
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distribution function F(a, %2,-+- , Xn) with corresponding characteristic function 
o( ,t,---, tn). Let G(x) be the distribution function of (a,X1 + --- + anX,)/ 
(b:X1 + +++ + bnXn) where a, , 2, °-* ,Qn,61,b2,-+-+ , bn are real numbers. If 
P{ Si b;X; < 0} = O, then 


G(x) + Gla — 0) =|]— =f {tay = bi 2x), . - ’ flan _ b,x) } dt. 
T 


Proor. Note that 
(~ , \ 
| aA; X ss { : P 
PY) x <2x> = P{rX(a; — bz) X; < 0}, 
wOi Ag ) 


and let R,(~) = P{Z(a; — bix)X; < &} and x.(t) be the corresponding char- 
acteristic function. Clearly R,(0) = G(x) and 
x2(t) - o{i(a = biz), Pee , t(a, -— bnx)}. 
On applying (2) to R,(é) and setting = 0, the required result follows at once. 
If (3) is applied in place of (2), with K = G, we obtain 
G(x) + G(x — 0) = H(O) + HO — 0) 
(4) 4 | { v(t) — o{t(a — bx), --- , ta, — baz} di 


mm t 


We shall consider (3) and (4) when n = 2 and 


(" by | 0 
ax by 7 0 1 


Two cases will be treated separately; first, when X, , X2 are independent, second, 
when X,, X2 may be correlated. 

If X; , X2 are independent, and F(a,, ©) = F(a), F(©, x2) = Fe(2xe), with 
corresponding characteristic functions ¢1(¢), ¢2(¢) then (1) becomes 
(5) G(x) + Gz - 0) =1- =f d1(Doo(— ta) iy. 

7 t 
while (4) becomes, taking H = F 
(6) Ga) + Ge -— 0) = 4 f b(t) — dda( =e) oy 
1 t 

Cramér [2, p. 46] proves, for X, , Xz independent and F.(0) = 0, that 


wv. . lL [* el) — o()d2(—t) 
G(x) = sal. can: mann dt 


under the following conditions: 


(i) X, and X; have finite means; 
i) fla < « 
1 
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If X; , X2 may be correlated, then (1) becomes 
(7) G(z) + Giz — 0) =1- +f oll, —t) oy. 
7 t 
while (4) becomes, taking H = F, 
(8) G2) + Ge — 0) = 14 (0 = oft, —te) 


Professor P. L. Hsu, in a course of lectures attended by the author at the 
Statistical Laboratory, University of California, gave the following result of 
Cramér, which was stated thus, using the above notation: 


1 ” pot) ea xed 3 
ot 


(9) G(x) = oni Le 


provided f _ xe | dy < oo, where F,(0) = 0 


and x,(t) is defined above expression (4). 

The following corollary is obtained from (1) according to well-known theorems 
concerning differentiation under the integral sign: 

CoroLLary. Suppose o(t1, f,-°-:, tn) is the characteristic function corre- 
sponding to X;, X2, +--+: , Xn, and G(z) ts the distribution function of 


(a,X, T+ + OnXn)/ (1X1 + -+- + baXn); 


(10) We) = A | a est dt, 
kel dt; tpent (a,—Bz2z) 


in every interval in which the integral converges uniformly. 


If n = 2, and 
aq by 1 0 
(* ") _ ( ). 
then 


vw. . 1 f Peel, t) 
(11) G(x) = af Ea 


Cramér (3, p. 317, exercise 6] states the following result: 


If F(u,%) = [ [. "f(u, v) dudv, and F,(0) = 0, then 


G'(2) = J f | ee} dt, 


if the integral is uniformly convergent with respect to x. 
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Geary [6] has shown that if F(x , x2) = [ | f(u, v) du dv, F.(0) = 0, and 


A(t, v) = | e"” f(u, v) du, then 
Rit ] [ 0¢o(t ’ By) 
we * "| ol Jena 


(i) ¢o(4, &) = Oforh = +~, 


1) | dy | yr(t, we dt = | at [ y(t, ye dy. 
0 = nn & 


provided 


Formula (1) can be employed in the case n = 2, X;, X2 are independent, and 


ay by 1 0 
a, de 7 0 | 


to obtain closed expressions for the distribution functions of ratios in which the 
variable in the numerator and that in the denominator may have any one of 
the following four distributions: Binomial, Rectangular, x’, Normal. In the 
case of the four ratios with the binomially distributed variable as the denomi- 
nator, a translation must be made to ensure positiveness of the denominator. 
For the four ratios with the normally distributed variable as denominator, the 
distribution function obtained is approximate; and the approximation is good 
if P{X2 < 0} is sufficiently small (cf. Geary [5]). 


5. Distribution of the ratio (a; X,+ --- +a,X,)/(b, Xi + --- + 6,X,), with 
denominator positive or negative. The following theorem will be proven: 

THEOREM 2. Let G(x) be the distribution function of (a,X1 + +--+ + a,X,)/ 
(bX, +--+ + bnXn) where ay, d2,°+: , Gn, bi, bo, +++, bn are real numbers. 
If P{ =? bX, = 0} = 0, then 


G(x) + G(x — 0) -1-14 


(12) 
{ta —b:x),--+,t(a,—b,2)} +¢@ {t(a, —biz),--- ,t(a, — b,x) di 
t ’ 

where 


Se.) = If ins f proct—tte a 


2bez,>0 


@ (i,t,-°+ ,t.) = [| eas / pe ae dF (x1 ,%2,°°* Xn). 


Lopz, <0 





INVERSION FORMULAS 


Proor. Let R, ¢) = P{=b.X; > 0} -P{ 2X. (a _ b,x) <é& | >b.X;. > 0} 
+ P{=b.X; < O}-P{2Xi(a —brx> — é | DhiX- < 0}. 


Then R,(o) = 1, R,(— ©) = 0, and R,(é) is non-decreasing in £ and continuous 
on the right. Hence R,() is a distribution function (Cramér [2, p. 11]). It can 
be shown by a proof analogous to that used by Curtiss [4] that the characteristic 
function of R,(é) is 


@ {t(a, — byr), --- , tan — bnt)} + {tia — ay), --- , t(baw — an)} 
Since 2.(0) = G(x), application of (2) to R,(é) yields the required result. 


6. Inversion formulae for multidimensional distribution functions. The 
n-dimensional analogue of (2) will now be given, and will be applied to obtain 
inversion formulae for the joint distribution of several ratios. 

Let X1, X2,°-- , Xn have the joint distribution function F(x , %,-++- , tn) 
and the corresponding characteristic function ¢(1, 2,-°--, t&). Let 


ee se ***% ti.) 


be the characteristic function corresponding to the marginal joint distribution 
function of X;,, X;,,°°* ,X,j,, where the set ji, jo,°+: , jx iS @ permutation 
of k of the integers 1, 2,---,n. Note that 


O(t1, 2, -++ te) = d1,2,---, nl, bys: , te). 
The summation 7 FE, , &i, » «°° » &ni,), Which will appear below is 
in) 


(81, $2, **°, §n 


to be interpreted as follows: 
Defining bi; = &; if 1; = A 
-— a 0 if 1; = 0, 
then 7 F(éi, , +++ , €ni,) Will mean that the summation is to be taken 
(i1, t2, -*°, tn) 

over all binary numbers .2%2 «++ tn. 

Using the notation of the preceding paragraph, we can state the following 
theorem: 

THeoreM 3. Let Ao, Ai1,°°* , An Satisfy then + 1 equations 


n—r—l 
p> (” ; ") Aves = 1, A, = -1, (r = 0,1,2,---,n — 2), 
k= 


where ("") as usual, denotes the binomial coefficient. 


(—1)"*" ; Peta, ,*** > &ni,) = Ao > p f¢ af 


(F1,8a,°°°.8 kat (1t)* i, <igces+ <ik 


(13) 
dt; dt, --+ dty 


-exp{ —t(t +--+ titi) } Pirie in plas +> te) hile +++ le 
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Proor: Since the theorem is already proved for n = 1 (section 3), and since 


1 P iltrErt++++tatn) dt; dt. +++ dt, 
(w2)” ff f é o(t . te . ’ a) "-hiaek 


=f, I= f sgn (a — &) sgn (a — &) - 


* Sgn (%. = En) dF (xy a Zn), 


the theorem could be proven by induction. The result is obtained more quickly, 
however, by noting that it suffices to consider the case of independent 
re ie 

It may be remarked that if (£:, &,---, &:) is a continuity point of 
F(a, , 22, °** , tn), the left-hand member of (13) becomes 


(—1)""'2"F(é ’ ae . En), 
and also that differentiation of (13) yields 


ae. eee » En) 


0&1 0 2 = > 0€n 


= (3) tf a f e tekites +take) o(h ; lo Ah a £.) dt, dt. cee dt, , 
2r 


in every n-dimensional interval in which the integral converges uniformly. 
This agrees with well-known results concerning Fourier inversion formulae. 
An inversion formula for the joint distribution of p ratios 


ai; Xi + Ari Xe +: HF Oni Xn ; 
bi X; + bes Xo +. + Dni An ; 


(14) 


= 1,2,---,p (l< pn), 


can be obtained from (13) by a method similar to that applied in section 4. ‘The 
following theorem holds: 
THEOREM 4. Let 


} | Yan X 


GG, &,-°: ,&) = P 


Lajp X; 
S&;° 


i és i « 
sox x, ‘2, En? 


and $(t:, t,°°:, tn) be the characteristic function corresponding to X, 
Xo,°°+,Xn. Then, if P{2bx~nX; < 0} = 0 (kK = 1, 2,--- p), 
(—1)?"" 2 Gb, ’ ois ee Eni,) 


(i, 42° tp) 


=A g thay, — bay, &, 
+BdhvaBaff ~felbm-mo 


= | dt; dtp +++ dt, 
ee Og Be) ip etnies 
? » i(anj, in 8i)9 b& os. 
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The following corollary is a generalization of (10) and follows by differentia- 
tion of (15): 
CoroLLARY. Suppose G(x, %,°-: , tp) ts the joint distribution function of 
the p ratios 
ay Xi - $ aX, 
bi X; + ee + Dnj Aa , 
and o(t; , 2, -* + , tn) ts the characteristic function corresponding to X,, X2, °° 
then, if Pie bijX; < 0} = 0, j = 1,2,--- p- 


anG (Erbe «+ Ep) _ (2 ) { ins ¢ 
) Bake == Ae Bei) 


; a” b(t, te, +++ tn) 
{> Dia Dye cee br» arp tem E15 (anj—bests) dry dro 


in every p-dimenstonal interval in which the integral converges uniformly. 


7. Joint distribution of ratios of quadratic forms. Let X,, X2, - 
have the joint probability density function 
is (det B)} e -}rBr' 
(2r)”/2 


f(x) 


where x = (2,22. °:: ,2,) and B isa positive definite symmetric matrix. Sup- 
pose @ is a positive semi-definite symmetric matrix of rank r < n and 
lL,, l.,--: , Ly is a set of symmetric matrices. We wish to obtain the joint 
distribution function G(& , & , +--+ , &p) of the p ratios 


XL, X' XL, X’ XL, xX’ 
XOX’?  XQN?’'’ XOX? 
where X = (X,, X2, °°: , Xa). 

The existence of such an orthogonal matrix S that SQS’ = I ”“ , where J” igs 
the diagonal matrix having the first r diagonal elements equal to unity and the 
rest equal to zero, is well-known. Let X = YS, C = SBS’, M; = SIS’, and 
note that C and the V/; are symmetric matrices. Also 


: ae if ™ 
Gg ‘ & ge se Se) a | \yT y’ < £4 eo eee ‘yyw yp < e,} ’ 


where Y = (Y;, Y2,--- , Yn) has the probability density function 


_ det C)! sey: 
gly) = (Qn)"?2 € re 
and y = (y1,Y2,°°° > Yn). 

Suppose the L; mutually commute in pairs. Then so do the M; ; for M;M; = 
SL;S'SL,;S' = SLL,’ = SL;L;S' = SL;S'SL;S’ = MM; ; since S is orthog- 
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onal. Hence, there is an orthogonal matrix U which simultaneously reduces 
each M to diagonal form; that is, N = UMU’ is a diagonal matrix (cf. Wey] 
[8, p. 25)). 


Let Y = ZU, D = UCU’, so that 





FAT r72 
ZN;Z' = De mi Z5 
j=1 


/ 2 y r72 
a ji 4; Mi» 13 
cane = a ee ee 
G(é , £2 ’ ? é,) P sr x — S 6; , b Sey" x < | , 
where v$” = lify <1; 
= OQifj>r, 


and Z = (Z,, Z2, ---: , Zn) has the probability density function 
(det D)} : 4zD2’ 





h() = 


(2r )nie : ? 
where z = (2, Z2,°°* , 2n)- 
We can now apply the results of section 6. If (h, &,---, tn) is the character- 
istic function corresponding to the joint distribution function of Zj, Z3,--- ,Z* 


it is clear that 
; det D ; 
bin *** ply a Epes E , 
Wh sb EE: ~ ath | 


where T is the diagonal matrix whose diagonal elements are t;, 2, --- ,tn. Ap- 


plying (15), with ¢ = y, we obtain, since G is obviously a continuous distribution 
function 


(—1)?*" 2? Gli, &, °°: , &) 


an 4+ 2 Gp nei adh” fy 12 may — vi &), 


, lu, --- dw 
e nee @ > Wins, wa HP gy) $B = _ 
i=1 W, We 
Ax, 
4 + at f 
° > (wi)* i, a> <: 
det DP doids dv 





-* cain ienas k ’ 
ie (r) Wi We +++ We 
det \ dag — 215.8 wi(us;, — vg &,) | 
lanl 


where D = [das] and das is the Kronecker delta. 


It is, of course, evident that a result analogous to (17) could be obtained, by 
considering p ratios 





INVERSION FORMULAS 


XIyX’—XIaX’ Xp X’ 
mr’ Mr’”6|lC CCU Re 


where the 2p matrices L,, Le, --- , Lp, Qi, Qe, --+ , Qp are symmetric and 
mutually commute in pairs, and Q, , Q2., --- , Qp are positive semi-definite. 

In the case p = 1 in (17) and for special classes of matrices L; , Qi, B the cal- 
culus of residues may be employed to obtain closed expressions for the distribu- 
tion of 


XIy X’ 
XQ, X’° 


Formula (17) can be applied to obtain the joint distribution of serial correla- 
tion coefficients with different lags. The author plans to incorporate these 
results with those mentioned at the end of section 4 in a forthcoming paper, 
written jointly with Roy B. Leipnik. 

The author wishes to acknowledge the valuable criticism of Professor H. Lewy, 
and especially the constructive advice and suggestions of Roy B. Leipnik. 
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THE FACTORIAL APPROACH TO THE WEIGHING PROBLEM! 


By O. KeMpPpTrHORNE 
Iowa State College 


1. Summary. The weighing problem is discussed from the point of view of 
factorial experimentation. The paper contains a brief description of the frac. 
tional replication of the 2" factorial system. It is shown that optimum designs 
for the weighing problem may easily be obtained with this approach. This 
approach is valuable in indicating the structure of weighing problem designs, and 
the limited conditions under which such designs can give results of value. 


2. Introduction. Considerable attention has been given recently to the 
problem of weighing a number of light objects on a scale [1, 2, 3,6]. The problem 
was originally proposed by Yates in his paper on complex experiments [4] as an 
example of a factorial experiment in which interactions between the factors tested 
would not be expected to exist: that is, the weight of say two objects could be 
assumed to be the sum of the weights of the objects weighed separately, after 
taking account of any necessary zero corrections. Such a situation is compara- 
tively rare in biological research when, for example, the effect on yield of a parti- 
cular crop from the joint application of two nutrients is usually different from the 
sum of the effects of separate applications. In recent years attention has been 
given to the use of fractional replication in factorial experiments [7, 8, 9] and it 
is proposed in this paper to consider the weighing problem from this point of 
view. 


3. The 2” factorial system. A full description of the 2” factorial system was 
given by Yates in his technical communication The Design and Analysis of 
Factorial Experiments [5]. Yates was particularly concerned with the analysis 
of such experiments and with the evolution of systems of confounding in order 
to reduce the number of plots in each block. The following brief account is 
given in order to facilitate the discussion of the weighing problem. 

In a single replication of the 2" system all combinations of n factors each at 
two levels are tested. With three factors, a, b, c, for example, the following 
eight combinations are tested: (1) a, b, ab, c, ac, be, and abc, where (1) denotes the 
control, a the application of treatment a only, ab the application of treatments 
a and b, and so on. A set of seven independent comparisons between the eight 
test results is given formally by the expansion of the formula 


Hat1lb+lce+ 1D), 





1 Journal Paper no. J 1548 of the Iowa Agricultural Experiment Station, Ames, Lowa. 
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where at least one of the signs is negative. If, for instance the first sign only is 
taken to be negative, a formal expansion gives the expression 


tiabe — be + ab —b+ac—c+a-— ll}, 


and this contrast of the observations gives the effect of the factor a averaged 
over the presence and absence of the factors b and c, which is denoted by effect A. 
Similarly taking the negative sign in the second bracket only, we get the average 
effect B, ; 
B = labe — ac + ab —a+be —c+b-— ll}. 
Taking negative signs in the first and second brackets we obtain the interaction 
AB 
AB = habe + c + ab + (1) — ac — be — a — Bb}, and so on. 
The definition of effects and interactions may be presented very simply in 
geometrical terminology, by representing the treatment combinations as points 
of an n-dimensional lattice, each axis of the lattice having two points at unit 
distance apart. The control treatment will have coordinates (0,0,0,--:,---, 
0), the treatment consisting of a only will have co-ordinates (1, 0, 0,--- , 0) 
andsoon. The effect A is then the difference of the mean yield of the treatments 
corresponding to the points lying on the hyperplane 
y= 0, 
and the mean yield of those represented by points lying on the hyperplane 
m= Be 


The interaction of two factors a and b, represented by the axes x; and 22 respec- 
tively, will be obtained from the difference of the mean yields of those plots for 
which 


X + rt = 0, or 2% + a = 2, 
and those for which 
m1+am=1. 


The extensions to the above for three-factor and higher order interactions are 
simple. The interaction of factors a, b, and c, which are represented by coordin- 
ate axes 21, 22, and 23, is given by the difference between the mean of plots 
represented by points for which 


X1 + Xe + 23 = Oor 2, 
and those represented by points for which 
t+ % + % = lord; 


in other words, it is the difference of the mean yields of those plots for which 
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ti + 2X2 + x3 = O (mod 2) 


and of those for which 


41 + X2 + x3 = 1 (mod 2). 





Each effect or interaction is then defined as the mean difference of two sets of 
plots, each set being represented by points on parallel hyperplanes, and the planes 
of one set of parallel hyperplanes lying between the planes of the other set. It is 
necessary to specify only the direction cosines of the hyperplanes in order to 
specify the effect or interaction, and the usual terminology for effects and interac. 
tions follows, in that the interaction of factors a, b, c, for example may be repre- 
sented by the symbol ABC. 

In the same way as effects and interactions are defined in terms of the yields 
of the several treatment combinations, the expected yield from each treatment 
combination may be expressed in terms of the mean level of yield and the true 
effects and interactions. If the full set of combinations of the factorial scheme 
is tested, the best estimate of each true effect and interaction is the same func- 
tion of the observed yields that the true effect or interaction is of the true yields, 
This fact is one of the advantages which follow from the use of the full factorial 
scheme. 

We are not concerned here with factorial experiments in which the factors have 
more than two levels, but when the number of levels of each factor is the same 
prime number, effects and interactions may be represented by products of powers 
of the symbols for the factors. In the case of two factors (a, b) at three levels, 
for example, the main effects may be represented by A, B, and the interactions 
by AB and AB’, each symbol referring to the two independent contrasts between 
three sets, each of three plots. 

As an example of the use of the above representation, we may consider con- 
founding, that is, the arrangement of the treatment combinations in blocks in 
order to reduce the experimental error. Suitable arrangements are such that 
contrasts between the blocks represent high order interactions which the experi- 
menter is confident will be of negligible size. 

If treatment combinations for which 


0X + ae%e + +++ + Qn tn = VU (mod 2) 
and for which 
B12, + Box +--+ + Br tn = O (mod 2) 
are arranged in a particular block, then the coordinates of the treatment combina- 
tions in this block also lie on the hyperplane 
(a1 + Bi)ti + (a2 + Be)te + +++ + (Qn + Bn)tn = VU (mod 2), 


where the coefficients (a; + 8) must be reduced modulo two. If, therefore, the 
treatments are arranged in blocks so that two comparisons are block contrasts, 
then the generalised interaction of these contrasts is also a block contrast. 
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4, Fractional replication. The principle of fractional replication follows 
very simply. Suppose only those treatment combinations whose yields all occur 
either in the positive or the negative part of a particular contrast are represented 
in the experiment, that is only those combinations represented by the points of 
the lattice for which say 


Oy 2% + ete + +++ + antn = O (mod 2). 
Then the comparison between the yields of those plots represented by 
B12) + Bete + --- + Brita = O (mod 2) 
and by 
Bi 21 + Bote + ++: + Bax, = 1 (mod 2) 
will be identical with the comparison between the yields of plots represented 
by 
(a; + f1)%1 + (a2 + Bolte + --- + (an + Bn)Xn = O (mod 2) 
and by 
(a1 + Bi)t1 + (a2 + Be)t2 + +++ + (an + Bn)tn = 1 (mod 2). 


The former of these two comparisons may be represented by the symbol 
dich? --- 2°", and the latter by af'**tag?t® ... aa"t?= where m1,--+, 2a 
are no longer coordinates but symbols for the n factors, which satisfy the relations, 
2? = 1, if a = 0 (mod 2). The equivalence of the two comparisons may be 
obtained by the use of an identity relationship in the symbols for the factors 


I = af'x3? +--+ 25" 


where J is interpreted as unity, and only those combinations whose coordinates 
(a, 22, °** » 2n) satisfy one of the equations 


ay 21 + aet%e + °° + antn = 0, or = 1 (mod 2), 


are represented in the experiment. If this identity relationship is multiplied 
by the symbol 24'x5? --- 2%" by ordinary commutative algebra, reducing the 
powers modulo 2 where necessary, we obtain 


+ + 
ah apf? re ah = 7{% Br) (aa Ba) 26. ap (tn tBn) 


It is more convenient to revert to the common use of capital letters A, B, C, etc. 
for effects corresponding to small letters a, b, c, etc. for the factors tested. An 
experiment in half-replicate is then represented formally by an equation of the 
type 

I = A*BC’ --- 


In such an experiment on n factors only 2” treatment combinations will be 
tested. Of the 2” — 1 independent comparisons in a fully replicated experi- 
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ment, information on one comparison is lost completely since only those treat- 
ments which appear in the comparison with the same sign are represented: 
the remaining 2" — 2 independent comparisons of a fully replicated experiment 
are identical in pairs giving 2" " — 1 independent comparisons. Each com- 
parison is then said to have two aliases and measures the sum (or difference, 
depending on which half of the treatment combinations are used) of two effects, 
an effect and an interaction, or two interactions. 

A quarter-replicated expériment can by the same process be represented by an 
identity relationship of the form 
I = A™BaI0N a a2 BB2c12 wee ae Artes) Gi t82) vi t+72) , 


It is useful in the evolution of fractional designs to note that the elements in the 
identity relationship form an Abelian group. 

Fractionally replicated experiments are formally identical with confounded 
experiments in that block differences may be regarded as additional factors in the 
confounded experiment. A 2” experiment arranged in 2” blocks, for example, 
may be regarded as a | in 2” design of a 2”*” experiment. Considerable care 
needs to be exercised in the use of fractionally replicated designs, but they have 
been found to be very useful in agricultural and biological research. 














5. The weighing problem. The problem of weighing a number of objects 
may be regarded as the problem of the estimation of the effects of a number of 
factors which do not interact. To take a simple case, consider the estimation of 
the effects of factors a, b, and c for which one complete replicate would consist 
of the combinations 

(1) a, b, ab, c, ac, bc, and abe. 


Suppose a half replicate design is used, based on the identity relationship 
I = ABC. 


The combinations tested would then consist of either the set {a, b, c, abc} or the 
set {(1), ab, ac, bc}. If the former set were chosen, the comparison estimating 
the effect A could also be ascribed to the interaction BC, that estimating effect 
B also to the interaction AB, and that estimating effect C to the interaction AC, 
as can be observed by multiplying the identity relationship by A, B, and C in 
turn. If the experimenter is confident that the two-factor interactions are 
negligible, then any effect given by each comparison would be ascribed to the 
main effect. 




















6. Discussion of a particular case. We give the derivation of a design for 
weighing a particular number of objects, say ten. Let the objects be denoted 
by a, b, c, d, e, f, g, h, k, 1. Then the total number of combinations which could 
be tested is 2", that is 1024, but as we are confident that interactions are negligi- 
ble, it is necessary only to estimate main effects. 

A fractionally replicated design must consist of a number 2” of combinations 
and this will be a 1 in 2°” design. A suitable fractionally replicated design 
consisting of 16 combinations will exist if it is possible to evolve an identity 
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relationship for a 1 in 64 design, such that each term in the relationship involves 
at least three letters. A possible identity relationship for such a design contains 
the numbers of the Abelian group obtained from all combinations of the elements 
1, ABC, CDE, EFG, GHK, ADL, and AFH, with the rule that the square of each 
letter is to be equated to unity. Each possible comparison may then be due to 
any of the 64 effects or interactions which may be derived from this identity 
relationship. In other words, each comparison has 64 aliases: in the case of ten 
of the comparisons, only one of the aliases is a main effect, and for the remaining 
five comparisons the aliases are all interactions of at least two factors. The 
actual design may be written down by finding combinations of the letters which 
have the same number of letters in common with the unit element and the six 
three-factor interactions. These are themselves a group consisting of all com- 
binations of unity and four combinations of letters. The sixteen combinations 
with an even number of letters in common with all the members of the identity 
group are found to be the following: 


(1) abdef,: acefl,* bedl, : 
abfgkl, degkl,: beegk, ° acdfgk, 
Sgh,° abdegh, aceghl, bedefghl, 
abhkl,’ defhkl, « beefhk,+ acdhk + 


The estimation of effects from the results of the sixteen weighings is particu- 
larly easy; the weight of object a will be one-eighth of the difference between 
those weighings containing a and those not containing a. There are ten such 
contrasts which estimate the effects, and the remaining five contrasts may be 
used to obtain an estimate of the experimental error. If o” is the variance of 
each weighing, the variance of the weight of a, that is, the effect A will be (1/8 + 
1/8)o” = (1/4)o”. The precision can be increased fourfold in the weighing prob- 
lem with a chemical balance by interpreting the absence of each letter as the 
placing of the object in the left hand pan and the presence as the placing of the 
object in the right hand pan. Each effect will then measure twice the weight 
of the corresponding object and the estimated weight of each object will have a 
variance of o’/16, that is, the same precision as if each had been weighed by itself 
eight times in each pan, or sixteen times in all. 


7. General case. ‘The rules by which fractional designs may be constructed 
have been exemplified above and the procedure is simple, though laborious in the 
case of a large number of objects. It does not, therefore, seem worth while to 
enumerate particular designs for the weighing of particular numbers of objects. 
A general procedure in considering the design for a particular problem is as 
follows. Taking the case of a number 7 of objects, the experimenter should 
form a rough idea of the order of magnitude of the experimental error, ¢ say, and 
decide what accuracy he requires for his estimates of the weights, a standard 
error say of s. Then if he weighs 2” combinations of the objects, the standard 
error of the estimate of each weight will be 2-””c in the case of the chemical 
balance. This serves to determine 2”” and therefore p, and it is then necessary 





















244 0. KEMPTHORNE 





to design a 2” experiment of fraction 2". Alternatively, a design of higher 
fraction which can provide estimates may be replicated the corresponding number 
of times. In the case of the spring balance the corresponding standard error ig 
2°?» 2. necessitating a design of higher fraction. 

Designs of the type described above have some useful properties: 

(1) the design automatically takes care of any bias in the balance, 

(2) the effects or weights may be computed easily as indicated above, 

(3) the effects or weights are uncorrelated, 

(4) all the effects are measured with the same precision, and 

(5) an estimate of the experimental error which is independent of the effects 

may be computed from the results. 

In considering the use for a particular problem of a design like the one discussed, 
it is important to understand completely the structure of the design. Such 
designs may well have real value for the weighing problem, but it is not easy to 
visualize other problems for which they would not give results capable of various 
interpretations. The use of the above designs depends on an assumption that 
interactions between pairs of factors are negligible, and this is generally not the 
case, for example, in biological research work, in which the experimenter may well 
be confident that interactions between three or more factors are negligible. In 
the particular case described in detail, there are only fifteen independent com- 
parisons between the sixteen results which will be obtained, and it follows from 
the identity relationship that the comparison giving the effect .A, also measures 
the two factor interactions BC, DL, and FH. If therefore the factor a has no 
effect and there is an interaction between factors b and c or the other two pairs 
of factors, the experimenter will conclude that the factor a has an effect. It is 
clear that. under these circumstances the experimental results are worthless. 


8. Efficiency of designs. In [2], Mood states for optimum designs such that 
when N weighings are made, the variance of the estimates of the weights are of 
the order of Ti in the case of the chemical balance and se in the case of the 
spring balance, where o’ is the variance of a single weighing. As indicated above, 
when JN is a power of 2, the fractional factorial designs result in the same 
variances. These designs are similar to those proposed by Kishen [6]. 

Mood dealt with the case in which the number of weighings was restricted, 
for example to be equal to the number of objects, and defined a best design as that 
which gave the smallest confidence region in the p-dimensional space for the 
estimates of the p-weights. 

To take an example for the weighings of 3 objects with a spring balance with 
no bias he suggests the following design: 


1 10 
X={|1 01 


01 1 
where the rows of the array refer to weighing operations and the columns refer 
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to objects. If the results of the weighings are y, , ye, and y; respectively, the 
estimates of the weights 6; , bz , and b; are given by the equation 


bi 1/2 1/2 —-1/2 Yi 
b}= | 1/2 -1/2 1/2] ly 


If oc’ is the variance of a single weighing, then the variance of each estimate will 
be [(1/2)? + (1/2)? + (—1/2)"Jo? = (8/4)o’: or if N (= 0(mod 3)) weighings are 
made by replicating the above system N/3 times, the variance of each estimate 
will be 90°/4N. The covariance of any two estimates is (—1/4)o” so that the 
square of the correlation between any two estimates is —1/9. The fractional 
factorial design will yield estimates which have a somewhat higher variance, 
namely 4c°/N. This increase in precision obtained in Mood’s design has been 
obtained at the expense of obtaining correlated estimatés which in addition are 
subject to any bias which the measuring instrument may have. It is doubted 
whether the use of such designs for any practical problem can be justified. 

It is of interest to note that the concept of fractional replication may be ex- 
tended to give designs requiring a number of weighings other than a power or 
two. For the weighing of 3 objects for example, a factorial design of fraction 
3/4 could be used: it could consist of a half-replicate based on the identity J = 
ABC, and a quarter replicate based on the identity 


I=A = BC = ABC. 


There is, however, little point in discussing such designs for the weighing problem, 
as their efficiency as measured by the total number of weighings needed to achieve 
a particular degree of accuracy is lower than for the designs described in this 


paper. 
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MULTIPLE SAMPLING FOR VARIABLES 


By Jack SILBER 
Roosevelt College 


Summary. A multiple (sequential) sampling scheme designed to test certain 
hypothesis is discussed with the following assumption: X is a random variable 
with density function P(x) which is piecewise continuous and differentiable at 
its points of continuity. Formulae are derived for the probability of acceptance 
and rejection of the hypothesis and for the expected number of samples necessary 
for reaching a decision. These formulae are found to depend on the solution of 
a Fredholm Integral equation. Explicit solutions to the problem are obtained 
for the case when P(x) is rectangular by reducing the fundamental integral 
equation to a set of differential-difference equations. Several examples are 
given. 























1. Introduction. A multiple sampling scheme is here proposed which is 
based on cumulative sums of random variables. Bartky [1] has developed a 
theory of multiple sampling for attributes when the attribute can take only two 
values with probability (p) and (1 — p) respectively. Formulae are there de- 
rived for the probabilities of acceptance and rejection of the null hypothesis and 
for the expected amount of sampling necessary for reaching a decision. In 
this paper the same type of formulae are developed for the case of variable samp- 
ling where the underlying probability law for the variable is given by a piecewise 
continuous function for which derivatives exist at its points of continuity. 

The whole theory of multiple sampling is closely related to Wald’s [2] theory 
of sequential tests. ‘The fundamental difference is that in the latter, probabil- 
ities of errors of the first and second kinds are assigned, and acceptance and 
rejection criteria derived therefrom, while in the former the problem is solved in 
reverse order. There the acceptance and rejection criteria are assigned, and 
probabilities of eventual acceptance and rejection derived. For different param- 
meter values, these are the probabilities of making errors of the first and second 
kinds. 

The problem presented here is similar to that given by Wald [3] in his paper 
on cumulative sums. In the present paper we waive the restriction that the 
expected number of items necessary for termination of the cumulating process 
be given explicitly as an integer. Since the theory given here is from the point 
of view of multiple sampling, the language appropriate to that theory will be 
used. 


2. The sampling scheme. Let X be a random variable with probability 
density function P(x) which is piecewise continuous. One variate, say 1, 
is selected and if z, > 6b, the hypothesis (for example the null hypothesis with 
respect to the mean) is accepted, and if 2; < a, the hypothesis is rejected. If, 
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however, @ < x, < b, another variate z2 is selected. In the latter case similar 
criteria with respect to x; + 22 determine whether the hypothesis is to be accepted 
or this method of sampling continued. Or more formally, let 


S. = 2a (r = 1, 2,3, ---), 
where the cumulative sums S, are formed sequentially as follows: for any integer 
r the cumulating process is terminated by acceptance of the hypothesis if S, > b 
and rejection if S, < a, but, ifa < S, < b another variate 7,4, is selected and the 
sum S,4: formed. The acceptance and rejection criteria are then applied as 
above. No attempt is made here to indicate the choice of the acceptance and 
rejection criteria. 


3. The probability of acceptance. If at the rth unit the hypothesis is neither 


accepted nor rejected, then it must be true thata < S,< b. Let us denote the 
probability that this condition holds by 


b 
(3.1) | Y,(S,) ds, ’ 


where Y,(S,) is the probability density function for S, in the above described 
sampling scheme. The probability density function for S,,; would then be given 
by 


b 
(3.2) , Y+410S8,41) _ | Y,(S,)P(Sy41 rac S,) dS, . 


The probabilities of accepting or rejecting the hypothesis on the rth trial are 
respectively 


+o a 
(3.3) [ vusjas,, [/ ¥.ts) as,, 
b —o 


and therefore the probabilities for eventual acceptance or rejection are given 
by 


(3.4) P= 2X [ ysas,, Pr=d | Y.(S) ds. 
fal © 


fl 


The probability that a < S, < b cannot exceed the probability that a < T, = 
a + te + x3 + +--+ + 2, < bona single sample of n variates, that is Pr (a < 
S, <b) < Pr(a < 7, < b). For distributions with positive variance, it can 
be shown that the right member of the above inequality approaches zero as n — 0. 
Therefore, the process of sampling as outlined above will eventually lead to 
acceptance or rejection of the hypothesis. See Wald (3, p. 284] for a direct proof 
that the probability that the left member of the above inequality holds for 
n = 1, 2, 3, --- is zero. 
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Consider the linear integral (Fredholm) equation 
b 
(3.5) V(x) = Yalz) + [ Pe — s)¥(s) ds, 


where Y;(x) = P(x) and assume a solution of the form 
(3.6) Y(x) = Yi(x) + AY2(x) + WY3(z) + ---. 


That solutions, in power series in \, of the Fredholm equation exist when the 
kernel P(a — s) and the function Y,(z) have finite discontinuities is well known 
and the theory has been expounded by several authors. (For example see 
Goursat [4].) If the power series in A is substituted in the integral equation we 
obtain 


Y,(x) + A¥2(zx) a Y;(z) +--+. 


Y,(x) +a [Yi(s) + Y¥2(s) +d” Yi(s) + -+-|P@ — 8) ds 
(3.7) 


b b 
rie +% | Y,(s)P(2 — s) ds +» | Tibi ~ ade 


b 
+x [ VsG)P@ — 9) dete. 
Equating coefficients of like powers of \ we see that 
b 
(3.8) V(c) = [| ¥,a()P@ — 8) ds, (r = 2,3,-++). 


This, however, is the probability distribution for S,, r = 2, 3, --- under our 
sampling scheme, and therefore from the equations, 


2 b 

(3.9) Y(x) = DANY, (x) = Vila) + | P(x — s)Y(s) ds, 
r=1 a 

we have that the probability of eventual acceptance for A = 1, is 


20 b 20 

(3.10) > [ vtsy as, = [ ¥@) ar. 

Thus our problem of finding a formula for the probability of eventual acceptance 
or rejection of the statistical hypothesis under the above sampling scheme re- 
duces to that of finding a solution of a linear integral equation. 

The argument in this section has, of course, been entirely formal. However 
from the general theory of integral equations we know that the series solution 
(3.6) converges uniformly for \ < 1/M(b — a) where P(x) < M, since P(z) 
is a probability density function. In equations (3.4) and (3.10) we give solutions 
for \ = 1 and of course we assume that M(b — a) < 1. Since (3.6) is uniformly 
convergent the interchanges of integration and summation in (3.10) and (4.3) 
in the following section are allowable. 
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4, The expected amount of sampling. Since 
b 
(4.1) | Y,1(S,-1) dS,-1 


is the probability that the rth sample will be reached, then the probability that 
on the rth sample, the hypothesis will be either accepted or rejected becomes 


b b 
(42) [ Ya») ds. - [ ¥.(8,) as, 


that is, the first term in this expression gives the probability that no terminating 
decision is made on the (r — 1)st sample and the second term gives the proba- 
bility that a like decision is made on the rth sample. The difference of the 
two then gives the probability that a terminating decision (acceptance or rejec- 
tion) will be made on the rth sample. The expected number of units sampled 
will therefore be 


b oo b b 
E 1— I P(x) dx + z r I Y,-1(S,_1) dS,-1 ie | Y,(S,) as, | 
a r=2 a a 


«© b b 
«Ze / Y,(S,) dS; = 1+ / T yiwe 


r=l1 r=1 


b 
1 + | Y(x) dz. 


Thus, the amount of sampling expected before a terminating decision is reached 
also depends upon the solution of the integral equation. We proceed to discuss 
the problem when P(x) is given by a rectangular distribution. 


5. Reduction to differential equations when P(x) is rectangular. Consider 
the integral equation 


b 
Y*(2) = P*(z) + / P*(2 — )Y*(2) dl, 


a 


PX(2) = 2, —cSz—-a< +e; 
(5.2) 2c 


= 0), z—-a>c or s~-a<—¢, 
and in the integral equation 
ata-c<z<b+a-e. 


The parameter a is restricted to the values — c < a < c for the following reasons. 
The rejection criterion a cannot be greater than c + a for, if so, rejection will be 
automatic on the first sample. Similarly the acceptance criterion b must be 
greater than —c + a for otherwise, acceptance would be automatic on the first 
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trial. If a > c then, rejection can never take place if it does not take place on 
the first trial for in this case all z > 0. Similarly, if a < —c then, acceptance 
can never take place if it does not take place on the first trial for in this case gl] 
z <0. Furthermore, in obtaining solutions of the integral equation, we wil] 
take a to be > 0. This inequality is no real restriction since solutions for nega. 
tive a can be obtained by symmetry from the solutions for positive a. 


If we let x = z — a then 





b 
(5.3) ¥*(c +a) = Px +a) +r] Pee +a—d¥*i a, 


(5.4) Y(x) = P(r) +X [ P(x — &Y*(d) dt, 











where 
P(x) = a, 
(5.5) 2c 
ey 
Now let s = ¢ — a, then 
b—a 


Par) +f eee 


—o 


Y (a) 


Il 


(5.6) 


b—a 


= Pe) +f Pe — a — YO) ds 


—a 


We have thus transformed our equation to one in which P(x) becomes sym- 
metrical with respect to the line x = 0. Furthermore, the probability of accept- 
ance becomes 


(5.7) P,= | Y(x) dz, 
b—a 
and the expected amount of sampling becomes 
b—a 
(5.8) E=1+ / Y(x) dex. 


Also, x now has the following range:a — c < «4 <b-+c. If we now make the 
transformation x — a — s = y, then 


(5.9) Y(~) = P(x) + [. P(y)Y(a — a — y) dy, 


and the following cases present themselves. 


Ife —a< —corx —b> +2, then Y(x) = P(z), since P(y) = 0, 
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Ife-b<-—-c<24—a < +e, then 
(510) ¥(z) =P@)+2[ Ye-a- dy, 


where 
a-cx<xrcsate when b — a > 2e¢, 


(5.11) 
a-c<z<b-—e when b — a < 2c, 


Ifz-—b < —c < +c¢ < x — a, then 
+e 

(6.12) Ye) =P@) +2 [ Ye-a-yay, 
where 
(5.18) a+te<2x<b—candb —c > 2. 
If-c<2—-b<2-—a< +e, then 
(5.14) Vic) = Pic) +> f v@-a-wady, 

2c z—b 


where 
(5.15) b—-c<2<a+candb —a < 2e. 
If -c <zx—b < +c < x — a, then 


. er 
(5.16) Y(x) = P(x) + = I - Y(x — a — y) dy, 


where 

b—c<2<b+c when b — a > 2c, 
(5.17) 

a+e<2<b+cwhenb — a < 2c. 


Transforming back to the variable s, we have for the case b — a > 2c, 


z—a+c 

Y(x) = P(x) + » | Y(s)ds for a—-c<x<a+te, 
x z—a-+e 

P(x) + > | Y(s)ds for a+e<a<b-e, 


b—a 
= P(r) +2 | Yim te b«06e e04e 
and for the case b — a < 2c, 


z—a+e 
Y(x) = P(x) + > | Y(s)ds for a—c<x<b-—e, 


b—a 
(5.19) = P(r) + 2 | Y(s)ds for b—ce<2ax<a+te, 


b—a 
= Pe) +2 | Y(s)ds for ate<2<bte. 
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In all of the above equations, the integral is a continuous function of z, a, a, b, ¢ 
while P(x) has a discontinuity at x = +c and x = —c, the jump at these points 
being of amount 1/2c. The function Y(x) will therefore be such that 


¥(—c + 0) — ¥(—c — 0) = 1/2, 
(5.20) 
¥(c — 0) — Y(e + 0) = 1/2c. 


If we now differentiate the above sets of integral equations with respect to x we 
obtain the following sets of differential-difference equations for the case \ = |, 
If b — a > 2c, 


1 


Y'(a) = — Yiw@—ate) for a-—cSzSate 


2c 
1 
2c 


(5.21) {Yaa—a+ec)—Y(@—a-—c)} for atcecxrcbdb—¢ 


2s 
2c 
and, ifb-—a< 


Y(a — a — c) for b—-c<a<bt+e 


Y'(x) : for a-c<2< 
for b—c<2 


=-5 Y@-a-0) for atesxrcb+b, 


the derivatives holding for all points except at x = —c and xz = +c. 

Although a technique has been devised to solve the above equations for finite 
a and b, mathematical difficulties of a computational character are encountered 
when (b — a) is made large. Note that there are only three essential parameters 
in the above problem since c can be taken as the unit of measurement. In the 
technique illustrated by the following examples, a has been fixed as has (b — a), 
i.e. the solutions shown in the examples below are general only insofar as one 
parameter is concerned. The essential feature of the technique is that the range 
of Y(x) has been further subdivided so as to make its points of discontinuity end 
points of subdivisions of its range, and thus Y (x) becomes continuous from the 
right or left in every subinterval of its range. 


6. Example I: b — a = 2c,c = 1,a = 0. In this case x ranges from (a — 1) 
or (—c), whichever is smaller, to (b + 1) or (+c), whichever is larger. If 
—c <a -— 1, then Y(z) = P(x) for -e < x <a — l,andifb + 1 < +e, 
then Y(x) = P(x) forb + 1<a2< —c. For zx between a — 1 and b + 1 the 
domain of the differential-difference equations is divided as follows, where a 
is now restricted to the interval —1 < a S 0. 
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(6.1) Yi(z) = 4Y¥ige(a + 1) where for i = 1, a-l<2z< -l, 
t= 2, —-l1l<2z <a, 
+= 3, a<2z<0, 
a= 4, 0O<2t<a+tl; 

Y((xz) = — 4Y,_2(a — 1) where for 7 = 5, a+1l<2xr< +1, 
+ = 6, +l<2z<a+t 2, 
t= 7, at+2s2z2< +42, 
a= 8, +2<2<a+3. 


The above are the equations corresponding to (5.21) for the given example. 
Differentiating the equations for 7 = 3, 4, 5, 6 and making certain obvious 
substitutions we obtain the following second order differential equations, 


(6.2) Y% (2) = —+Y.(z), i= 3, 4, 5, 6, 


where the domains for x are as in (6.1). If we solve the equations (6.2) and sub- 
stitute in the remaining equations in (6.1) we obtain the following set of equa- 
tions, 


Y;(z) = A ius sin 4(x + 1) _ Buse cos (a + 2) + K;, i= i, 2, 
(6.3) Y.(x) = A; cos 3x + B; sin 3z, 4 = 3, 4, 5, 6, 
Y,(z) = —Aj_esin $(x — 1) + Bi_2 cos }#(x@ — 1) + Ki, 4 = 7,8, 


where again the domains are as in (6.1) 

From continuity considerations we have the boundary conditions 
Ya-—1)=Y¥sa+3)=0, Yi(—-1)-—3= ¥2A(-1), Y2(a) = Ys(a), 
Y:0) = ¥40), Yaa+1)= Ys@+ 1), Yos(1) = Ye(l) + 3, 
Y.(a + 2) = Y,(a + 2), Y,(2) = Y,(2). 

These boundary conditions yield certain relationships between the constants. 
The equations so determined, however, do not form a consistent set of linear 


equations in the A;, B;, K;--:-. If we integrate out the equations (5.18), 
sectionally, the following relationships between the constants are obtained. 


A; = Aizesin } — By42 cos 3, B; = Bj42 cos $ + Bie sin 3, t = 3, 4, 
Ky = —(A4 — Ag) sin 3(a + 1) — (Bs — By) cos $(a + 1) 
=~ $+ Bh — B+ Ki, 


K; = As — As + Ke, Kg = Ag sin 3(a + 2) — By cos $(a + 2), 
(6.4) 


Bz; = 4 + Be t+ Ky — Ag + As + (A, — As) sin (a + 1) 
+ (B; — B,) cos 3(a + 1), 
Ag = As + Kg + (Ag — As) sin 3(a + 1) + (Bs — By) cos $(a + 1), 


kK, = —As sin 5 + B,; cos 5. 








254 JACK SILBER 





From these equations it is easily seen that Ay = A; and K, = Kz = K; = Kg. 
Furthermore, the following set of consistent linear equations is obtained, after 
several simple manipulations and substitutions. 


(sin Ma +2) +sin 5 - sin i} Ag 
_ {00s a py {cos 4(a + 2) + sin 5 — COs sa = 4 
(6.5) {- sin 3(a + 2) + cos #(a + 2) — C08 5 — + gin i} Ag + “ “} B; 


+ {sin 4(a + 2) + cos (a + 2) — cos 5 Os cos 3} B, = 
) 








{cos 4} As — Bz; + {sin 3} Be = 0, 


All the other constants can be obtained from the solutions for A, , B; , Bs in 
(6.5). Letting A equal the determinant of coefficients in (6.5) and using the 
relationships (6.4) we obtain the following solutions: 






A = 2 — 2sin 3} — cos }, 

AA, = ${cos } — cos a/2 - sin (a + 1)} = AAs, 

AB, = 3{sin 4} — sin a/2 - sin }(a + 1) + cos 3} —1}, 

AA, = #{sin 1 — cos 3 + cos a/2 - cos #(a + 2)}, 
(6.6) AB, = }{sin $(a + 2) cos a/2 — sin $ — cos 1}, 

AB; = 3{1 — sin a/2 - sin }(a + 1) — sin 3}, 

AAs = }{cos } — sin’ }(a + 1)}, 

AB, = 3{sin } + sin 3(a + 1) - cos 3(a + 1) —1}, 


















Ak, = 1 {cos § sin 4(a + 2} = AK, = AK; = AKg. 
If we now integrate Y (x), equation (6.3) sectionally, i.e. from the left end point 
to the right end point of each sub-interval of its range and then add up appro- 
priate areas, we obtain the following formulae for the probabilities of acceptance 


and rejection and for the expected amount of sampling: 
{cos }(a + 1) + sin a/2 — cos a/2 + AKo}, 
— {2 — cos} — 2 sin} + sin 3(a + 1) 


— cos }(a + 1) — sin a/2 + AK3}, 


= F {cos a/2 — 2 sin a/2 — sin }(a + 1)}. 
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7. Example II: a = 1,c = 3,b —a=4. In this case, as in the previous one, 
Y(z) = P(x) for - 3 <x <a—3whena—c=a—3 < —3andifb+c= 
a+7 <3 then Y(z) = P(zt),a+7<2<3. Fora—3 <x <a+7where 
a takes on only integral values between —5 and 3, we have the following set of 
differential-difference equations: 


Yasi(t) =  &Yass4e(x + 2), j@ ~3, ~2, -1,6; 

(7.1) = 0, j = 1,2; 
= — $¥a4j;-u(x — 4), j = 3, 4, 5, 6. 

If we integrate the above equations for 7 = 1, 2, substitute in the equations for 


j= —1, 0, 5, 6, integrate, and then substitute in the remaining equations, we 
obtain the solutions 


Yori(t) = te Aatssa(z + 2)’ + Act i422 + Ac;, j = -—3, -3; 


= gAa+j42% + Aas; = —1,0; 


J 
= Aas; j 1, 2; 
J 


— pyAay;-2(a — 4)” — $Aayi4t + Acs;; = 3,4; 
= — tAajur + Aas;, j = 5,6. 


As in the previous example we now use (5.22). Integrating out (5.22) sec- 
tionally, certain relationships between the Aai;, 7 = —3, —2,--- 6, are ob- 
tained. These yield 


Aoi = gs {12P.-1 + 12P, + 39Po41 + 9P 42}, 
Aai2 = gs {12P.4 + 12P, + 11Pa41 + 37P 442}, 


Aga = — 5; (4Po1 + 4Po + 13Pan + 3Pais] 
+ y5{228P.1 + 60P. + 55Pay1 + 17P 42}, 
=— — {12Po1 + 12P. + 11Pas: + 37P ais} 


+ 7¢s{60P.1 + 228P. + 55Paq1 + 17Po42}, 
where P,,; is the value of P(x) fora +j<4sa+j+1,j = —3,---6. 
All of the other constants can be found in terms of those given in equations 


(7.3). If we now integrate (7.2) sectionally and perform several simple manipu- 
lations, we arrive at the sii formulas: 


he 3a — a 3 l 
Pe — x, Pay + + Agu + -—..- Aat2 + To fe = a 4° 


3a — 89 9a — 131 , 
> Pari 916 Aa + ~ 916 * A ate + prAa1 + PsAc, 


2a — il — 13 
— 1 + 12 Aa += 12 Aat+2 + Aa + As. 
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Although P.,;,j = —6, —5, —4, 7, 8, 9, have not appeared in previous deriys, 
tions in this example, they have been inserted in the above formulas to cover the 
cases in which a — ¢c > —corb+c<e. 

It should be mentioned that Kac [5] obtained the distribution of n (the ex. 
pected amount of sampling) by a process similar to that given in this paper. |; 
is also interesting to note that the present paper could have been written entirely 
in the language of problems in Random Walk. 

The author has also worked on the case in which the distribution P(z) is ty. 
angular and parabolic. In these, as in the case of the rectangular distribution 
discussed in this paper for b — a large, the equivalent differential-differeng 
equations are of large orders making the computation of solutions extremely 
tedious. As Kac [5] pointed out, the task of obtaining solutions in closed form 
for the case when P(x) is the normal law is extremely difficult. 
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ON THE CHARACTERISTIC FUNCTIONS OF THE DISTRIBUTIONS OF 
ESTIMATES OF VARIOUS DEVIATIONS IN SAMPLES 
FROM A NORMAL POPULATION 


By M. Kac 


Cornell University 


1. Summary. An explicit formula for the characteristic function of the 
deviation 


the, a > 0, 
N kel 


is derived for samples from a normal population. For a = 1 one can calculate 
the probability density function but the result does not seem to be in complete 
agreement with a recent formula of Goodwin [1]. 


2. Introduction. Let Xi, X2,--- , X, be independent, normally distributed 
random variables each having mean 0 and variance 1. 
Let, as usual, 
a Xi + Xe + ay +X, 


X — ee ——oe ee — ’ 


n 


and denote by Y,(a) the deviation 
, : 1 n ‘ a 
(1) Y.(e) = —- >), | X,— X| ’ 
N kal 
The purpose of this note is to show that 


F,(€) = Efexp (éY .(a))} 


l ” P —r a] n(Elz|*%+-nzr % 
= Facvear® [| [Leaner ae | an 


It is easy to check that for a = 2 one obtains the well known expression 


Qit re 
e-2 


Moreover, if a = 1 one can actually find the probability density of Y,(1). The 
resulting expression is fairly complicated and it strongly resembles an expres- 
sion recently obtained by Goodwin [1]. Except for the relatively simple case 
n = 3, it does not seem easy to verify that our formula is equivalent to that of 
Goodwin. 

Although deviations corresponding to values of a different from 1 and 2 are of 
little practical value the explicit formula (2) may be of some interest. It is 
also hoped that the method of derivation may prove useful in other cases. 
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(2) 
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3. The derivation of (2). We start with the observation that 
X and Y,(a) 
are statistically independent (see e.g. Daly [2]). 
Denote by 

E*{ |X| < «, exp (#€¥,(a))} 


the integral of exp (2Y,(a)) extended over that portion of the sample space in 
which | X | < «. Thus the conditional expectation E{exp (i€Y,(a)) | |X| <¢ 
is given by the formula 








Efexp (itX.(a)) || X| <q) = ZU sis ai ae . 





Because of the independence of X and Y,(a) we have 


(3) E{exp (#Y,(a))} = ah 7 { = 














For the sake of simplicity we assume now that a > 1 and note that 
exp (#Y x(a)) — exp (| Xe") SED || xe me FT 
o HANS mez po 
Thus, on the portion of the sample space where | X | < ¢, we have 
“exp (Yala) ~ exp (# Di xii) |< BE Ki to 
and consequently 


Et | X | < €, exp (7€Y n(ax))} er BY | X | < €, exp (@xix, ry 
< He { |x <6 D(x to}, 


Clearly Bt |X| < «> (|X| + \ , approaches 0, as ¢ approaches 0, 
1 
hence by (3) 


Bs | X | < cen (#> xsi") 


(4) E {exp (Y,(a))} = lim omen <8 ca 
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Using the fact that 
|X| <« 
* [22 exp (ink) dn = 4, |X| = 6 
Tle 7 7 
|X| ><«, 


we obtain easily 


a |X| < sexn(# |X: r)} 


a net ef exp Yo @| Kal" + oxo ban 


Tle 


5) 


The justification of interchanging of the order of integration (from - to «) and 
the operation E can be made quite simply (see e.g. Kac and Steinhaus [3]). 
Notice now that 


B{ exp © (& | X; |* +x} 


o lz [ exp (- 5 exp £ (| [* + nz) as | = ¢,(€, ) 


and that ¢,(é, 7) is absolutely integrable in (— ©, ©) as a function of 7. 
Thus, as e > 0, 


(6) : [ 7 = 7 on, 0) dn ~ = [. onlé, 0) dy. 


Furthermore (as « — 0) 


(7) Prob {|X| <e} +20 VE 


and combining this with (6), (5) and (4) we get 

sai 1 - 
(8) E{exp (#Y,(@))} = y= [ on(é, 0) dn. 
This, of course, is equivalent to (2). 


4. Density function of the mean deviation. If a = 1 one can obtain an ex- 
pression for the probability density f,(8) of Yn(a). We note first that 


C) 2 . 
[- exp (- 3) exp £(@ | 2 | + we) ae 


” 22 
=n I exp (- tS) exp i(€ + )x dx 


© 22 
+ nf exp (- “f) exp u(& — nx dx = n{R(E + n) + REE — n)} 





R(u) = [ exp (- = exp (tux) dz. 


Using (2) (with a = 1) we obtain 


7 ee. ae “Ts (n k n—k 
F,(€) = J nr/2n)"* [. P @ R'E + yn) R'E — » | dn. 
n. 


Let us first look at the summands corresponding to k = 0 and k = We have 


[ R'(E — ) dn = [ R"(n) dn = [_re + n) dn. 


Now, R(n) is the Fourier transform of 


eo 2£<¢ 


xp { — s> 
exp ( =). x 0, 


and hence R”(n) is the Fourier transform of the convolution 
Ceo ee ee = FM(Z), 
n 


It is easily seen (using integration by parts) that 


R(n) = O (4) 


for large | 7 | and hence for n > 2, R”(n) is absolutely integrable in (— 0, ~). 
It follows (classical inversion formula) that 


| R"(n\ dy = 2x¢™ (0). 


Since for n > 2, ¢(zx) is continuous and ¢ (x) = 0 for x < 0 we 
have ¢”(0) = 0. Thus 


grt 8 6s — - 
F,(é) = (/ae) Xu will RY(E + n)R™“(E — 2) dn. 
It is fairly easy to check that 


[ R*(E + 9)R"“(E — 9) dn = [ exp (itr) (5) rr (3 


so that 


ie mnt . s “., eS (n (k) [2% (n—k) [ @ . 
F,,(é) — (\/2n)"*? [. exp (iéx) 2 c )r (3) c (5) dx. 
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Finally, 


Fe (*) ww 4) (n—k) .) 
HB) = (gaye 2, (e) #” (5) s°* (5). 
[ have not been able, except for n = 3, to verify directly that this formula is 
identical with that of Goodwin. 
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A FUNCTIONAL EQUATION FOR WISHART’S DISTRIBUTION 
By G. Rascu 
State Serum Institute and University of Copenhagen 





1. Introduction. The sampling distribution of the moment matrix for ob. 
servations from a multivariate normal distribution was given by Wishart in 
1928 [1]. This proof involved rather advanced multidimensional geometry but 
since then two analytical proofs have been given: one by Wishart and Bartlett 
in cooperation with Ingham by the use of the characteristic function [2] anda 
second by Hsu by induction with regard to the dimension of the observa- 
tions, [3], 

In the following section is given a new derivation of the form of Wishart’s dis. 
tribution in which a fundamental property of the multivariate normal distribv- 
tion is utilized, viz. the invariance of the distribution type against a linear trans. 
formation. In section 3 the same principle is used for evaluation of the constant 
and determination of the moment matrix in the multidimensional normal 
distribution. 













2. Derivation of Wishart’s distribution. Let’ 
(1) x = (m%1,-°-- , 2%), 
denote a k-dimensional normal variate with the mean vector 0 and the distribu- 
tion matrix 
(2) 
viz. 
VAS) | -ixex* 
(3) p{x} = VAS). -iner, 
(/2n)* 
® is symmetrical and positive definite. 
Now consider n observations of x: x;, --* , Xn, which are stochastically in- 
dependent. Their joint distribution is 


(4) P{m%i, °°: Xa} = (MARY : eo Pex r Oxy 


The estimation of @ is based upon the moment sums 










= (vis), 












ni; = LLL; ’ 





1 Notations: Lower case latin and greek letters are scalars; boldface capital latin and 
greek letters denote matrices, and boldface lower case letters row vectors. * means trans- 
position. A (A) stands for the determinant of the square matrix A. 
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which form the symmetrical, positive definite matrix 
(5) M = (mj;) = 2x; x. 


In order to derive the distribution of M the straightforward procedure seems to 
be to transform the distribution of the sample (x: , --- , Xn) to a distribution of 
M and some other variables which then should be integrated away. As such, 
the transformation, 

(6) x, = u,M?, M2u?u, = 1, 


might serve. The matrix 
(7) 


contains 7/: elements linked together with * * ” relations; (U) symbolizes 


(1 ae > ') k of the elements taken as independent variables. 


For the purpose of introducing M in the exponential term in (4) we shall 
define the ‘“‘double dot multiplication” of two matrices: 


(8) A-- B= (ay) «+ (bi) = DO Do ashiy, 


(i) (7) 

for which we notice the rule 
(9) A-- (BCD) = C-- (B*AD*). 
As obviously 

xbx* = Yo; 7.7; = &-- (x*x), 
we have 
(10) >x,ox; = &-- M, 
and accordingly 


a(M, (U)) |’ 


Ae) cleo | a tes Xn) | 


(11) piM, (U)} = (co 


) , ; ; ' 
\) denotes the jacobian of the transformation. On integrating with 


‘ 
where — 


a ) 


respect to (U0) we obtain 
(12) piM} = (V/a(@))" - o**™ - o(M), 


where ¢(M) is independent of &. From this it follows that p{x:,--- ,Xn|M} 
is independent of ®, i.e. M is a sufficient statistic for ®. 

In order to determine the mathematical form of ¢(M) we shall apply an arbi- 
trary linear transformation to the original variates: 


(13) x, = xA. 
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The new variates x, are obviously normally distributed about 0 with the dis. 
tribution matrix 


(14) ® = AdA*. 













Therefore the distribution function of the new moment matrix, given by 


(15) M = A*M’A, 
is 
(16) p{M’} = (Va@’))” - *"™ o(M). 


On the other hand the transformation from M to M’ is a linear one, the jacobian 
of which therefore is a constant depending on A only: 


{ 
(17) aan = ¥(A), say. 
Consequently, t 
(18) p{M’} = VA) - «*™ - o(M) - | ¥(A) |. 
The two expressions for p{M’} must be identical, and as 
(19) A(®’) = A@)A*(A), 





and 


(20) 

















&’ -- M’ = (A®A*) -- M’ = (A*M’A) --@? = M-- , 
it follows that ¢(M) satisfies the functional equation 
(21) | ACA) | “e(M’) = ¢(M) - | ¥(A) |. 


Now, since the transformation M = (AB)* M’(AB) may be carried out in two 
steps, ¥(A) also satisfies a functional equation 


(22) ¥(AB) = ¥(A)y(B). 
Furthermore, if A is a diagonal matrix it is easily seen that 
(23) ¥(A) = (A(A))*™, 


and this relation holds generally. In fact, considering the case where the normal 
form of A is a diagonal matrix: 


A= TDT ", say, 
we get 
W(A) = V(T)y(D)¥(T”) 
= (A(D))*" (TT) 
= (A(A))", 


and by analytical continuation this is seen to be true for any A. 
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Now, inserting this result in the functional equation (21) and taking for A the 
real symmetrical square root of M so that? M’ = 1, we readily obtain the 
solution 


Ge) e(M) = (A@M'))"**-9(1). 

_ ‘it follows that 

(25) p{M} = y:(n)(A(@))”? - et ™ . (A(M))**?, 

where yi() = ¢(1) is a constant which may be determined in various ways (cf. 


for instance Cramér [4]). 


3. Other applications of the linear transformation. It may be noticed that 
the linear transformation also leads to simple derivations of two fundamental 
‘ properties of the normal multivariate distribution itself, viz. determination of 
the constant and the relation between the moment matrix and the distribution 
matrix. 

Let 


(26) p{x} = y@)-*™, 
and transform by 
(27) x = yA. 


The new variable obviously has the distribution matrix (14) and the constant 
y(®’). But on the other hand direct transformation of (26) leads to 


0(x) 
d(x’) | 
y(®) | A(A) | eo *™", 


Pix’) = y(@) 


and therefore we must have 

y(®’) = y@) | A(A) |. 
For A = &' we get ’ = 1 and consequently 

y(@) = VA) - r(1), 
where obviously 


1 


y(1) = (aa) 


Considering 


M() = [ x*xp0x) dx, 


2 Exists because M is positive definite: Let M = ODO* where O is orthogonal and D 
the diagonal form of M; then M+ = ODiO* is real and symmetrical. 
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the same transformation gives 


M(s) = | A*x* xAp{x'} dx’, 


= A*M(®’)A 
which for A = ®° leaves us with 
M(@) = @’)” 


because M(1) = lL. 
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THE DISTRIBUTION OF A DEFINITE QUADRATIC FORM 





By HERBERT ROBBINS 


University of North Carolina 


Let X1,--- , Xn be independent normal variates with zero means and unit 















variances, let a; , --- , dn be positive constants and define 
aSvts... 422 

(1) UO. 9 X 1 + + 9 X n> 

(2) F,(z) = Pr([U. <a], fa(x) = F(z). 

Setting 

(3) a= (a; eee ar 

and using the convolution formula we may show by induction that for z > 0, 
“ C) k 

(a fla) = ahah SY Om) 


tao T(3n + k)’ 


aha a) 
t=0 '(3n + k + 1)’ 


(5) F(z) 
where for k = 0,1, --- 


, 2)...1(4 a 
©  aeet 5 Mtb ThtD,, 


e > ‘ 
ippesdigmk =otal e+ tnlay'++-a,” 
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In particular, if a; = --- = a, = 2, then using the known distribution of x’ with 
n degrees of freedom we have 


gi! —jr gi gi! 2 


. Po (-~ é.(—2)* 
Jalz) = Tan) 2 T(4n) Ss ] po a3 2, T'(in + k)’ 
so that 
4 Tmt FP yp G+) TH +P 
2 kT (3n) PF bgtee-btgnh tl---ta! P 


and therefore 


Ta + 3)---TG, + x" T(an + k) 
6y-bo + +46 gunk “hh. eee ——s kIT ($n) 


(7) 


Now in the general case let 
(8) a = min {a,--- , Gn}; 
then from (6) and (7) we deduce that 


| 


(9) |_ex(—2)" | < 2/2)" 


iTGn + k)| > S TGnye’ 
with a similar inequality for the general term of (5). 

‘It is difficult to evaluate numerically the coefficient c, by a direct application 
of the definition (6). We shall therefore give a method by which the c, may be 
computed easily. We shall assume in what follows that the a; are distinct. 

Let Y;,---, Yn be another set of variates with the same joint distribution 
as the X; and independent of the X;, and set 


(10) Vn = SX $e + EK SV 4 + EYE, 


(11) Gon(x) = Pr [Von < x}, Gon(x) — in 
Then by the convolution formula, 
z CS) k k 
(at) = _ = ate 31> eens) 
(12) Jan(x) I fala y fly) dy a & oe | {z Ci cash T(k + n) F 
But we may show directly that, setting 


(13) qg: = II (a; — a)” 


ins 
we have 


(4) gala) = ("De gat tet = (— eS na a (=2)" 


k=O k! 


Equating coefficients in (12) and (14) we derive the fundamental formula 


(15) a CCei = a” > “—e (k =0,1,--- 
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Thus, defining 


(16) 2P, = a" i qi eo , 


t=1 


we may write 
k 
(17) > Cex; = 2P,. 


From (6) or (17) it follows that 
(18) Co = . 


Thus we may solve (17) ene for the c;, in terms of the P;: for7 = 0,1, --- 


c; 
(19) Co = Po; — { eves a Ht CeCe Fes Hpi + a} 


Costa = Pajyr — {€1C2j + C2C2ja + +++ + CjCj41}. 


When the n constants q:,--- , dn have been computed we may compute the 
P; by (16) and then the c; by (19) successively as far as desired. The inequality 
(9) gives an indication of the number of terms of the alternating series (4) or (5) 
which are necessary to secure a desired accuracy. A sharper result than (9) 
should certainly be possible when the a; are well separated. 

The foregoing method may be modified to cover cases in which some of the 
a; are equal, the formulas (16) being replaced by the appropriate limits as the 
a; approach equality. 

We shall now derive an expansion of f,(x) and F,,(x) in terms of x’ distribu- 
tions. Let us set for x > 0, 


n+k—1 ewe 


— = _ 4k x 
(20) fin(x) dX ( L) di. gintk T(an + k) 


or, equivalently, 


§n+k—1 —}a 


a a x é 
3! (5 *) = > es Qk T(dn + k) 
(21) 


a 7 1 e * oo (—1)": ‘ 1, (ae)" T'(3n) 


Qin i = r'(43n ‘+ - k)’ 
where the d; are to be determined. It follows, after some reduction, that 
(22) Qc) =a"a"'e"" > > d; dij \ ‘ (—z)" : 
k=0 | i=0 | ak T(4n + k) 


But we may write (14) in the form 


— ~awn o~t, —zla —2—k _— «| oes xy tt 
(23) genx) = a" | . aqeay "(a — ai)’ faHT(k + 1)' 
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Equating coefficients in (22) and (23) and setting 
‘ —(k+1) n P 
(24) 20), =a ys qi a; ™ (a — a;) iii Pe 
i=l 


we obtain the relations d) = 1 and 


k 


ya dy = 2Q« 


1=0 


from which the d; may be computed as in (19). Equation (20) or (21) then 
gives the expansion of f,(x) in a series of x’ frequency functions. The corre- 
sponding expansion of F,,(x) is then 

cd i gin te—t ete 


F(x) = (—1)' dy - 
(26) ” 2, yd Jo aih** Tn + k) ™ 


k--0 


F (5 r) > (-1)'d [ a 
aa. = i Segoe ag at. 
2 M.S Sg QP (Aan + k) 


By direct comparison of (4) and (20) we may establish the following relations 
among the c; and d,: 


h 1 — 
dy = (—1)* D (-a)’ ‘) cj, 
j=0 


k—j 


£. fi i in 
ch = a Zz v - ‘) d;. ° 


j=0 


(28) 


‘ ° 2 . ‘. 
These may be used if both the power series and the x° series are desired. 
From (6) we see directly that 


lic _ 
(29) a == Da;’, 
i=l 


_ 


and from (28) it follows that 
(30) d, = 4 anb,, 


where we have set 


(31) by = 15 2. az \ ae (a, BB, ey ; 


Since by a well known inequality b; > 0 it follows that d, > 0, with equality only 
if allthe a; are equal. If we denote by h,(x) the frequency function of 
4a(X} + --- + X%) then 


in—1 —iz 


a a a é 
¥ sh (5) = dart: 
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and hence (21) becomes 
(33) 


which is significant for x near 0. 
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EXACT LOWER MOMENTS OF ORDER STATISTICS IN SMALL SAMPLES 
FROM A NORMAL DISTRIBUTON 


By Howarp L. JoNEs 


Illinois Bell Telephone Company 


1. Summary. Exact means in samples of size < 3, and exaet second moments 
and product-moments in samples of size < 4, are given in Table 1 in terms of 
a for order statistics selected from the normal distribution N(0, 1). The deriva- 
tion employs multiple integration and some general properties of the moments, 


TABLE 1 
Expected values of lower moments of order statistics, x; > Zizi , 
in samples of size n from the normal distribution N(O, 1). 


Moment | n=3 


En) =z 8/(@V/x) 
E[x2] 0 
Elz) 1+ — 
E|\z3] 1— 3/7 
E|zx,22] V/3/(2r) 
E[x,23] —V/3/r 
eel (2/3 — 3) 
E|x2xs] 2V/3 — 3)/a 
ot 1 — (9 — 2v/3)/(4z) | 
1 — V/3/n 

012 \/3/(27) 

013 ; (9 — 44/3) /(4x) 





2. Introduction. The usefulness of the lower moments of order statistics for 
determining the moments of the range and for other purposes is well established. 
In small samples, however, computation of the moments by quadrature is labori- 
ous [1]. The values shown in Table 1 should therefore be helpful in problems 
requiring the use of these moments for samples of size < 4, since the constant t 
has been evaluated to several hundred decimal places. Some of the methods 
used to obtain these results may also be useful in approximating or verifying 
the moments in larger samples. 
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3. Multiple integration. Let » random selections from the normal distribu- 
tion N(0, 1) be arranged in order of size so that 


aA em S *** SM. 
For samples of size 2, the means and product-moment are easily obtained from 
the general formula 
- / ata 
1 ze 


Elxiz*] = wf ff 
f(ar)f (ae) +++ flan) dai dxz +--+ dz, 


Zn- 


1 2 
f(x) = va a, 

E [zi] being the special case where h = 0. Multiple integration can also be 
used to find any product-moment, E [2;x;,:], for samples of size 3, the order of 
integration being changed at any stage where necessary. 

For the means in samples of size 3 and the product-moments in samples of 
size 4, the integrals reduce to double integrals which can be evaluated from the 
equation 


ra —(a2t24+b2¢2) Tv 
€ mv?’ dt dg = —. 
a J, te Sab 


This equation follows from the fact that 


“ £* ab —(a%¢2+b2¢2) 
—e mv"? dt, dte 
— 30 bte/a 


T 


is equivalent to 


1 1 
[ J ene, 
0 “pe 


while the function 


te 


has the symmetrical property that ¢ (4) = —¢(—&), whence 
| $(t2) dt, = 0. 


4. Some properties of the moments. The most obvious property of the 
moments of order statistics in samples from the normal distribution N(0, 1) is 
their symmetry; thus: 


E [x —E [xn-i+1), 
E [zi] E [2tn-:41); 


E [xixj) = E [2n-i41%n—j+1). 
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When sample values from any parent distribution are numbered in order of 
random selection, x; and x;,4; are statistically independent of each other, and 
the expected value of a product 2jz’; is the product of the expected values of a 
and x’. Numbering in order of size has the effect of increasing some expected 
values and decreasing others, leaving the sum of expected values of a given type 
unchanged, so that in general, 


n—1l n 
> dL (Elxiz')) = (5) E{xdlE[xol 
i=l j=i+1 


=t+ 







where 2% is a random selection. In particular, this equation holds for the special 
cases (k = 1,h = 1), (kK = 1,h = 0), and (k = 2, h = 0); so that in samples 
from the normal distribution N(0, 1), 






Gace d n 
> (Elea)) = in(n — (Ela) = 0, 
t=1 j=i+l 





- (E[x,]) = nElaxol = 0, 
i=l 





>» (E[x‘]) = nk\x5] = n. 
i=l 





The foregoing relationships lead immediately to the evaluation of EF [2,2x2| and 
E {x}] in samples of size 2. (The generalization of these relationships was sug- 
gested by Professor John H. Smith, whose unpublished manuscript on sampling 
from a rectangular distribution has also been instructive.) 

In samples from a normal distribution, the covariance of every order statistic 
with the sample mean is the same as the variance of the sample mean. This 
implies that the variance of the sample mean < the variance of any order sta- 
tistic, the ratio of one standard deviation to the other being equal to the co- 
efficient of correlation between the sample mean and the order statistic. To 
derive these properties, consider the linear function 


m = WyX1 + WeXe +--+ + WaXn 
of the order statistics X,, X2,--- , X, in a sample selected from the normal 
distribution N(u, o) with unknown yz and oc. Let 
ui = (X; — p)/o, $= 1,2,---,n. 
The conditions w; + we + --+ + wa = Land wa_iq, = w; are sufficient to make 
m an unbiased estimate of » with variance o” E [(witi + wate + +++ + wate). 
The w’s that make this variance minimum must satisfy the equations obtained 
by replacing w; with w,_-i4: , for > 3(m + 1), in the expression 
E[(wiar + wees + +++ + Wada)’ | + Mwi + we + +++ + wa — 1) 


and then setting the partial derivative with respect to each w equal to zero. 
This leads to 



























Zz Ww; Elz; xj] + 2. Ww; Elan x3] + A= 0, 1 “ 1 < n, 
jul j=1 







pecial 
mples 


| and 
5 sug- 
ipling 


utistic 

This 
T sta- 
1€ CO- 


To 
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where the summations include the terms E[z{] and E[x,-:4:'], respectively. But 
it is known [2] that the sample mean is the regular unbiased estimate of » with 
minimum variance. Setting each w equal to 1/n and combining equivalent 
terms yields 


. E[x;x,;) + 3nd = 0, += 1,2,---,n. 
i 
Summing from 7 = 1 toz = n, and employing the relationships discussed in the 
preceding paragraph, we obtain 
n+ 1nd = (), 
whence 
= —2/n, 


and 
ae Elx;x;] = 1, i= a 2, en, 
j=1 


where the summation includes the term E[z;]. This equation leads to the prop- 
erties mentioned at the beginning of this paragraph. The same equation can 
be used to evaluate E[xj] and E[z3] in samples of size 3 or 4 from the distribution 
N(0, 1), after the product-moments have been found. 
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NOTE ON AN ASYMPTOTIC EXPANSION OF THE nTH DIFFERENCE 
OF ZERO 
By L. C. Hsu 


National Tsing- Hua University, Peiping, China 


This note gives an asymptotic expansion of the nth difference of zero. It is 
known that the Stirling number S,,,, of the second kind is defined by 


(1) n'S,, = A'O' = PB (—1)"* (") x’. 


z--0 


We shall first show that the Stirling number S,,,,. can be expanded in the 
form 


2* ke! 


2k " oh " 
. én? E 4A), lb) Sl) ow) |, (t<k) 


n n n' 
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where fi, fe, --- , f: are polynomials in k and whose coefficients can be found 
by means of the following lemmas. 


The first lemma is due to B. F. Kimball, [1, (5.3)]. 
Lemma 1. (Kimball) Let ¢ be a real number such that n + q > 0, and let 
f(c) = x”**. Then we can write A"f(x) in the form 


(3) A*f(z) = fw + 4n) 1 += £. '). es + ~s W(m, n)|, 


m=1 





where the value of W(m, n) is given by 
(4) W(m,n) = Bom (—3n)/(3n)””, 






By"(x) being a so-called Bernoulli polynomial of negative order which was first 
defined by Nérlund [2). 


LemMMA 2. Let the sum of au(? *) products of k different numbers taken from 
the set (1,2, --- ,n) be denoted by S;(n). Then we can express it in the form 


(5) Sim) =D (-r,0(2 1), 
where the coefficients d4(k), \2(k), --- satisfy the recurrence relation 
(6) (kK + p)dAp-i(k) + p-d,(k) = A, (k + 1) 





with X%» = 0, Ar = 1 and Agai(k) = 0. 











Proor. Clearly, among all ) produets of (k + 1) numbers out of 


k ‘ 1 
(1,2, --- ,~), there are exactly (” k ') products containing the greatest factor n. 


The sum of these products is therefore n- S.(n — 1). Repeating this reasoning, 
we get 


(7) Srai(n) = n-Si(n — 1) + (n — 1)-Si(mn — 2) + --> + (K+ 1)-S2(K). 


Evidently, (5) is true fork = 1. Suppose now that it istruefork = k. Then 
the right-hand side of (7) can be written as 


: i ici > ie prry(” + ; rte = ') 


-> phon] os na ut ted 2) _ (. bet i) 








SC kt+p+1 i : n+p 
= (HD + apalt) + p- lly ~*~ i) 


The lemma thus follows by induction on k. 
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The number S;() may be called a Stirling number of the first kind. By the 
lemma just proved, it is easy to find 


s(n) =-3("F7)_("t1) 

sin = (2) (O24) +02 
w sin 0(*$4)— (°F) 4 (°$)-C21 
°) 1260(” + *) + 490 (” . *) 

9 8 
ae n+2 n+1 
so("77)+("5'). 

We shall see that in order to compute the coefficients of f,(k), fo(k), «++ , it is 


sufficient to compute the values of W(m, 7), \1(m), A2o(m), --- , (m = 1,2, ---, 2). 
Let f(z) = x"**. Then by lemma 1, we have 


a! Sank = | oH os | {1 a : ‘is : wom, n) | 


From the definition of S;,(n) it is easily seen that 
(nt k)\(n+k—1)---@+1) =n" + n™"S(k) + +--+ + 2Sia(k) + Si(k) 


Hence we may write 






ot < 


S:(n) = ae 









1 i(k) Salk) S.(k) 
4) Fie +3 |= JF, (+2 + +: +50), 


It is clear from Kimball’s paper [1] that 


. (wom, n)=)>> io W(m, n) + O(n”). 


m=1 n=l] 






Substituting, we obtain 


Santt = oe eal + = i. ‘) W(m,n) + on | 
{i + : = | 
sift + (ok,) Wom) + 06m) | 


™m=al 


{i f= bs + °) on) __ +0(n) |. 


mal p= \M + p/ (—1)°(—n)™ 











(9) 








ae 
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The last expression shows that the asymptotic expansion (2) can be obtained 
by computing the numbers A,(m), W(m, n) with 1 << p<m<t. For example, 
consider the case t = 3 and notice that [1, (2.13)] 


= nm LT i n nr n i n 
B n aati _ nine B n sees ee ae : R n agar 
' ( *) 12 ( *) 48 120 - ( *) 


and that A: = 1, A2(2) = 3, Ae(3) = 10, A3(3) = 15. Then by a straightforward 
calculation of the right-hand side of (9) and by comparison with (2), we find 


fuk) = 32K + b) 
(10) fo(k) = se (4k* — k’ — 3k) 


: (40k° — 60k’ — 2k* — 63k*° + 133k? — 48k). 


fx(k) = 810 


Finally, combining (2) with the well-known Stirling’s formula [3] 


te afoenl(*V 11424.) — 3 sr, 
~~ ae vie Onl + ton > aaen2 ~ sistone + 2M ) |p 


and noting (1), we obtain 


2\k n : ; ‘ 
cai) arom = ¥2en() (2) [1 +) 4 AO) 4 OO) + 00-4] 


k! 2 


where i(k), ge(k), gs(k) are polynomials in k, viz. 


1 - 
k) = — (8k° + 4k 1). 
qul 19 | + 4th + 


gk) = |. (e4k — 0k + 1). 
8 
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l 6 5 3 
as(k) = ~.—.. 2560k° — 3840k° + 832k' — 4032k 
51840 
+ 8392k° — 3732k — 139). 
The asymptotic formula of A"0"** just derived is much better than a result 
previously obtained [4]. Moreover, it may be noted that the asymptotic ex- 
pansion of S,,.4; may be made as sharp as desired, since in fact, for any pre- 
scribed ¢ > 1, A,(m) and Bz, (—4n), (lL < m < 2), may be easily computed by 
(6) and Kimball’s [1, (2.12)] respectively. 
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AN INEQUALITY FOR KURTOSIS 


By Louis GuTTMAN 


Cornell University 


1. Summary. It is well known that, if the fourth moment about the mean 
of a frequency distribution equals the square of the variance, then the frequencies 
are piled up at exactly two points, namely, the two points that are one standard 
deviation away from the mean. In this paper is developed a general inequality 
which describes the piling up of frequency around these two points for the case 
where the fourth moment exceeds the square of the variance. In a sense, it is 
shown how “U-shaped” a distribution must be according to its second and fourth 
moments. 


2. An inequality. Jet x be a random variable whose distribution has the 
following moments : 


(1) uw = E(x); 0° = E(x — yw); (ae + Lo = E(x — yp). 


a’ is non-negative for any distribution, and its positive square root will be denoted 
by a. Let 


(2) t= (x — p)/o. 
It will be shown that, if \ is an arbitrary positive number, then 
(3) Prob {1 — Ae Sf S1+ra} >1— rv”. 


If \ is chosen so as to make the left member in the braces positive, then ¢° is 
bounded away from zero, and (3) becomes: 


(4) Prob {+/1 — ha S|t] S$ V1 + ra} > 1-2”, (Aa < 1). 


For example, if a = .5 and A = 1/2, then (4) shows that the probability is 
greater than .50 that ¢ is either between .54 and 1.30, or between — 1.30 and —.54. 
Ifa = .2and\ = 3, then (4) shows that the probability is greater than .88 that 
tis either between .63 and 1.27, or between —1.27 and —.63. In general, the 
smaller a is, the greater the probability that ¢ is in a small interval around +1 or 
—1. In particular, if a = 0, then \ may be taken arbitrarily large, so that (4) 
shows that the probability is unity that ¢ = +1; this is the well known case 
referred to above. 
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3. Derivation. Inequality (3) is a special case of‘a slightly more genera) 
inequality which follows very simply from that of Tchebychef. Consider the 
function  — 1 + c, where c is an arbitrary real number. By using (1) and (2), 
it is seen that 


(5) B@-1+cP=a'+e°. 

Then, according to Tchebychef’s inequality, if \ is an arbitrary positive number, 
(6) Prob {( — 1 +c)? S$ Ne? + c*)} > 1 — Yd, 

or, 

(7) Prob {1—c-AVeFt+e es(sl—c+raAVeit+e} >1—-r% 


This is the general inequality that was to be shown. 

Inequality (3) is obtained by setting c = 0 in (7). 

Another special case is obtained by determining c so as to maximize the left 
member in the braces of (7). By differentiation, the maximizing value is found 
to bec = —a/+/d? — 1, for which (7) becomes: 


(8) Prob {1 — o8 $f $1 + a(6’ + 2)/8} > 1 — 1/6 + 0), 


where 8 is used instead of the notation 1/2 — 1, and denotes any positive num- 
ver. For the same probability on the right, (8) has the advantage over (38) of 
having 1 — af greater than 1 — Xa, so that the former may be positive even 
though the latter is negative. Inequality (8) starts the positive interval for ¢ as 
close to +1 as possible. On the hand, (3) provides the minimum size interval 
for ¢ from among all values of c that make the left member in the braces of (7) 
positive. 

If it is desired to have the positive interval for ¢ end as close to +1 as possible, 
then the right member in the braces of (7) is to be minimized. By differentia- 
tion, the minimizing value is found to be ¢ = a/+/)? — 1, and the minimum in- 
equality is: 


(9) Prob {1 — a(@’ + 2)/8 s #@ S$ 1+ a8} > 1—- 1/(6 + 1). 


4. Distribution Around yu. If the left member in the braces of (7) is negative, 
then instead of giving information about the piling up of probability of ¢ around 
+1 and —1, (7) provides a statement about the probability in an interval around 
u. Alternatively, this may be regarded as a confidence interval for wu. The 
minimum interval is given by (9); actually, it holds regardless of the value of the 
left member in the braces; another way of stating it is: 


(10) Prob {-V/1 + a8 StS V1 + ae} >1-— 1/@'+ D). 
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TABLE FOR ESTIMATING THE GOODNESS OF FIT OF EMPIRICAL 
DISTRIBUTIONS 


By N. Smirnov 


1. Editorial Note. The table presented on pp. 280--281 was originally pub- 
lished in [1]. It gives values of 


L(z) = 1 _ ° bm faa" gr ne (Qn)? 2 : or 
- vorl 


which is also derived in [2]. 

Let (X1,-°°* , Xn) be a sample of independent variables with the same con- 
tinuous cumulative distribution function F(x), and let N(z) be the number of 
X, which are < z. Fy empirical distributicn is meant the step-function 
F*(z) = N(z)/n. The maximum D, of the difference | Fa(z) — F(z)| is a 
random variable and L(z) is the limiting cumulative distribution function of 
nD, . If Dmn is the maximum of the difference | Fx(z) — F2*(z) | between 
the empirical distributions of two independent samples of sizes m and n, respec- 
tively, then L(z) is also the limiting cumulative distribution function of 
(mn/(m + n))"*Dmn . 

REFERENCES 
(1) N. Smirnov, ‘On the estimation of the discrepancy between empirical curves of dis- 
tribution for two independent samples,’’ Bulletin Mathématique de l’Université 
de Moscou, Vol. 2 (1939), fasc. 2. 


[2] W. Fever, ‘“‘On the Kolmogorov-Smirnov limit theorems for empirical distributions,”’ 
Annals of Math. Stat., Vol. 19 (1948), pp. 177-189. 
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TABLE of L(z) 
z L( 
.28 .000 
.29 .000 
.30 .000 
ol .000 
.o2 .000 
.30 .000 
.o4 .000 
.30 .000 
.36 .000 
oF .000 
.38 .001 
.39 .0O1 
.40 .002 
41 .003 
.42 .005 
43 .007 
44 .009 
.45 .012 
.46 .016 
47 .020 
.48 .024 
49 .030 
.50 .036 
51 .042 
.52 .050 
.53 .058 
04 .067 
.55 O77 
.506 .087 
57 .098 
.58 .110 
.59 .122 
.60 .135 
.149 
.62 . 163 
.63 177 
64 . 192 
.65 .207 
.66 .223 
67 .239 
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TABLE of L(z) 


Continued 

z L(z) 

69 .272 189 
.70 .288 765 
71 .305 471 
.72 .322 265 
.73 .339 113 
74 .355 981 

75 .372 833 
.76 .389 640 
77 406 372 
.78 .423 002 
.79 4389 505 
.80 .455 857 
81 A472 041 
.82 .488 030 
.83 .503 808 
.84 .519 366 
.85 .534 682 
.86 .549 744 
.87 .564 546 
.88 .579 070 
.89 .593 316 
.90 .607 270 
91 .620 928 
.92 .634 286 
.93 .647 338 
.94 .660 082 
.95 .672 516 
.96 .684 636 
.97 .696 444 
.98 .707 940 
.99 719 126 

1.00 .730 000 
1.01 .740 566 
1.02 .750 826 
1.03 .760 780 
1.04 .770 484 
1.05 .779 794 
1.06 .788 860 
1.07 .797 636 
1.08 .806 128 


.3o0 
.36 
.OF 
.38 
.39 
.40 
41 
.42 
.43 
44 
45 
.46 
47 
.48 





Z 


09 
10 
ll 
12 
13 
14 
15 
16 
17 
18 
19 
.20 
21 
22 
.23 
24 
25 
26 
27 
.28 
.29 
30 
31 
32 
oo 
34 
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TABLE of L(z)— 
| Continued 


L(z) 


.814 
.822 
.829 
837 
.844 
851 

.858 
.864 
.870 
.876 
.882 
.887 
.893 
.898 
.902 
-907 
.912 
.916 
.920 
.924 
.928 
.931 
.935 
.938 
941 
.944 
947 
-950 
.953 
.955 
.958 
.960 
.962 
.964 
.966 
-968 
.970 
971 
.973 
.974 


342 
282 
950 
356 
502 
394 
038 
442 
612 
548 
258 
750 
030 
104 
972 
648 
132 
432 
556 
505 
288 
908 
370 
682 
848 
872 
756 
512 
142 
650 
040 
318 
486 
552 
516 
382 
158 
846 
448 
970 
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TABLE of L@)— TABLE of L@)— —_—|| TABLE of L@)— 
Continued Continued | Concluded 
2 | L(z) | 2 
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.998 421 
-998 536 
.998 644 
.998 744 
.998 837 
-998 924 
.999 004 
.999 079 
-999 149 
.999 213 
.999 273 
.999 329 
.999 380 
.999 428 
.999 474 
.999 516 
-999 552 
.999 588 
.999 620 
.999 650 
.999 680 
.999 705 
.999 728 
.999 750 
999 77 

.999 790 
.999 806 
-999 822 
.999 838 
.999 852 
999 

.999 

-999 

.999 

.999 

.999 

.999 

.999 

-999 

.999 
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BOOK REVIEW 


Fundamentals of Statistics Truman Lee Kelley. Harvard University Press, 
1947; pp. xvi, 755. $10.00. 


REVIEWED By A. M. Moop 


Iowa Slate College 


First, a brief look at the contents: introductory matter, broad classifications 
of types of data, quantitative and qualitative aspects of data, construction of 
tables, charts, and graphs—200 pages; location and scale parameters, and 
moments—75 pages; normal distribution—30 pages; exact sampling distributions 
based on normal theory—5 pages; binomial distribution, goodness of fit tests, 
contingency tables, normal approximation to the distribution of the variance 
ratio, properties of Chi-square—20 pages; correlation and regression—150 pages. 

These first 480 pages constitute the essential part of the book and the part that 
will be commented on here. But there are 270 more pages, the content of which 
we shall merely note without comment. There is a chapter of 90 pages entitled 
‘Sundry Statistical Issues and Procedures”’ which discusses fi‘teen issues such as 
periodicity, time series, curve fitting, variance error of a coefficient corrected for 
attenuation, machine extraction of square roots, and sequential analysis. There 
follows a chapter of 40 pages devoted to no less than twenty-three topics in 
mathematics, topics such as: matrices and determinants, the square root trans- 
formation, expanding a table, spaces of three or more dimensions, and Fourier 
series. The remaining 140 pages contain numerical tables, references, various 
indexes, and a test designed to measure the adequacy of students’ mathematical 
preparation. 

This then is another book which deals with the descriptive aspects of statistics. 
Despite its title, it omits discussion of distribution theory, sampling theory, 
the theory of estimation, tests of hypotheses, or the theory of probability. The 
phrase ‘‘confidence interval”’ appears not once, I believe, in the entire 750 pages. 
The discussion of Student’s distribution is brief enough to be quoted in its 
entirety (page 284): ‘The t-distribution, shown through the courtesy of Dr. 
Philip J. Rulon, in Chart VIII II, is appropriate for interpreting the significance 
of means, differences of means, and of regression coefficients, for small samples— 
say N less than 15. It is the distribution of these statistics computed from small 
samples drawn from a parent normal distribution.” 

Thus the author denies any value to the developments in the fundamentals of 
statistics during the past twenty-five or thirty years. He does this not merely 
by implication but in so many words; referring to modern statistical inference, 
he writes (page 13): ‘‘A still greater weakness is that it is essentially a deductive 
procedure and relatively sterile in suggesting new courses—in inspiring creative 
inferences. It is fundamentally a method of proof and not one of invention.” 

282 
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He is therefore fully aware of his extreme position, and takes great pains to 
justify it. His thesis is that the main purpose of statistics is to suggest new 
hypotheses to the scientist. In developing this thesis he writes (rage 15): “‘The 
physicist observes seemingly irregular changes in z as y changes. He repeats 
his experiment, controlling more and more of the conditions, and repeats again 
and again, and, if successful, he reaches a law at the end of his work. He has 
been using statistics.’”’ But his discussion avoids certain relevant questions. 
Why does the physicist repeat the experiment? Why did he perform it in the 
first place? Did he suspect before he collected any data that x and y might be 
related ? 

At any rate, the opinion of most present-day statisticians is that the primary 
role of statistics in scientific research is statistical inference. This opinion is 
certainly well-founded in my own experience. Here at Iowa State College the 
Statistical Laboratory is intimately implicated in the research programs of all 
departments—physical, biological, and social. These scientists perform their 
experiments with a specific purpose in mind—usually the estimation of some 
parameters, sometimes the testing of a hypothesis. They never seem to seek 
in a collection of data some new hypothesis by artful selection between the mean, 
the mode, the geometric mean, the harmonic mean, and the median. 

It must be reported that, even as a book on descriptive statistics, it leaves 
much to be desired. The errors usually found in such books are to be found here 
as well as many more. There is the long discussion of skewness and kurtosis 
based on the false notion that moments are determined by the nature of the 
distribution in the neighborhood of the mean. Certain properties of the normal 
distribution are imputed to all distributions. Erroneous criteria for selecting 
amongst the many means are given. The universality of the normal distribution 
seems exaggerated; thus, for example, referring to deviations from regressions 
(page 364): ‘‘Since the quantities (%» — Zo) are ‘errors’ we may regularly assume 
them to be normally distributed.” Population parameters and their estimates 
are confused. The book contains a great many statements (like the final one in 
the section on the Student distribution quoted above) which are so carelessly 
written that they have to be counted as errors. Several of the derivations and 
arguments are also carelessly constructed; an extreme example of this appears on 
page 206: ‘Is the mean an unbiased statistic? M = (a1 + t+ 2% + °°: + 
tn)/N. Since the various 2’s are independent, there are just NV degrees of freedom 
and M is unbiased.” 

Students will likely have difficulty with this book. There is an air of arti- 
ficiality because of the omission of any discussion of population distributions and 
the notion of random sampling. Without any background of this kind it is 
hard to motivate the presentation, and the various topics become isolated. 
Moments are defined in terms of sample observations, and population moments 
are defined merely as the limits of these moments as the sample size becomes 
infinite. To introduce the mean, the author writes essentially: let us consider 
the function f(b) = [Zx°/N]'*. There is no pointing to the middle of a distribu- 
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tion function, or even a sample, or a histogram. The variance is introduced the 
same way; one considers the function = | 2; — x; \"/(N? — N). Technical terms 
are used without definition; for example, in the passage about the mean quoted 
above, the student suddenly encounters the word “unbiased” without definition 
or previous discussion and must infer its meaning from the context. 

Perhaps the best part of the book are three chapters on correlation and regres- 
sion. The idea of correlation is here introduced with the discussion of a numerical 
example, and several other topics are discussed in terms of examples. This part 
of the book is very exhaustive; every sort of correlation coefficient is discussed 
as is every sort of correction to such coefficients. But still the writing is careless, 
and there is some confusion of ideas. The worst confusion occurs because the 
distinction between normal and intraclass correlation is never brought out; the 
discussion hops back and forth between the two ideas with no hint that they 
are not the same thing. This part of the book, too, is in the style of statistics of 
thirty years ago; the emphasis is on correlation coefficients rather than regression 
coefficients. 





NEWS AND NOTICES 


Readers arc invited to submit to the Secretary of the Institute news items of interest 


Personal Items 


Dr. Leo A. Aroian of Hunter College has been promoted to an assistant 
professorship. 

Mr. Carl A. Bennett is now with the General Electric Co., Hanford Engineering 
Project, Richland, Washington, as an engineer in the Statistical Division. 

Dr. Arthur B. Brown has been promoted from an Assistant Professor to an 
Associate Professor of Mathematics at Queens College, Flushing, New York. 

Professor Maurice H. Belz has returned to the University of Melbourne, 
Carlton, Australia after having spent six months in the United States. 

Dr. Edward E. Cureton, member of Richardson, Bellows, Henry «& Co., Inc., 
industrial psychologists, is now at the United States Naval Air Station, Pensa- 
cola, Florida working on a project with the Navy. The object of this project is 
to improve ground school training, especially instructor training, in the Naval 
Air Training Command. 

Mr. Eric F. Gardner has accepted an assistant professorship at the School of 
Education, Syracuse University, Syracuse, New York. 

Mr. Lee S. Gunlogson, formerly with the Lumbermens Mutual Casualty Co. 
at Chicago, is now with the Marketing Services Division, Carrier Corporation, 
Syracuse, New York. 

Dr. Theodore E. Harris has accepted a position with the Douglas Aircraft Co., 
Santa Monica, California. 

Dr. Manuel O. Hizon, a former graduate student in the Mathematics Depart- 
ment, University of Michigan, is now with the Bureau of Banking, Manila, 
Philippines as Actuary-Examiner. 

Mr. Julius Lieblein, formerly in the Treasury Department, Washington, D. C., 
has transferred to the Statistical Engineering Laboratory, National Bureau of 
Standards, where he is working on problems in acceptance sampling and process 
control. 

Mr. Jack Moshman, formerly a tutor of mathematics at Queens College, 
Flushing, New York, has been appointed to the staff of the Department of 
Mathematics, University of Tennessee. 

Dr. Horace W. Norton, formerly with the U.S. Weather Bureau, Washington, 
D. C. as meteorologist, is now at Oak Ridge, Tennessee. His position there is 
to study the application of statistics to reliability of weighings and analyses in 
connection with accountability for source and fissionable materials. 

Mr. Emil D. Schell of the Bureau of Labor Statistics has been appointed Chief 
of the Mathematics and Electronic Computor Branch in the Office of the Comp- 
troller, United States Air Forces. 

Miss Bernice Scherl, formerly with the Schenley Research Institute, Inc., New 
York, has accepted a position as Statistician, Shell Oil Co., New York. 
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Dr. Irving E. Segal, who has been an assistant at the Institute for Advanced 
Study at Princeton, New Jersey, has accepted an assistant professorship in the 
Mathematics Department, University of Chicago. 

Miss Rosedith Sitgreaves, assistant statistician in the United States Public 
Health Service, has returned to her position in Washington after doing advanced 
study at Columbia University. 

Dr. John E. Walsh, who received his doctor’s degree in mathematics from 
Princeton University last October, is now employed by Douglas Aircraft Co., 
Inc. of Santa Monica, California. 

Mr. Winfred P. Wilson, a former graduate student at the University of Michi- 
gan, has accepted an assistant professorship at the University of Houston, 
Houston, Texas in the Department of Mathematics. 


a 


Announcement has been received of a new journal, The British Journal of 
Psychology, Statistical Section, which is published by the Council of the British 
Psychological Society. The editors are Professor Sir Cyril Burt and Professor 
Godfrey Thomson. The first issue has been published and later issues will be 
published as material warrants. Subscriptions and inquiries should be sent to 
the University of London Press, Ltd., Warwick Square, London, E. C. 4. 


—— a 


Announcement of Navy Department Joint Board of U. S. Civil Service 
Examiners 


Implementing its scientific research and development program both geo- 
graphically and in new fields of endeavor, the Navy Department is currently 
expanding three comparatively new, permanent laboratories in California. 
Heretofore, the Navy Department’s scientific centers have been concentrated in 
the eastern and eastern seaboard areas. 

Two of the laboratories have been established as the logical outgrowth of 
programs carried on by universities during the war. The Naval Ordnance Test 
Station, China Lake (formerly Inyokern), California, 160 miles from Los Angeles, 
was originally an activity of the California Institut2 of Technology. Its present 
program involves research, development and test work with ordnance equipment 
and explosives. The Navy Electronics Laboratory at San Diego, California is 
the outgrowth of work done by the University of California. It is concerned 
with research, testing and development of electronic control devices, detection 
equipment, instrumentation equipment and training aids. The Naval Air 
Missile Test Center at Point Mugu on the coast of California, 60 miles north of 
Los Angeles, was established when the need for an installation became apparent 
as the result of the Navy Department’s activities on guided missiles. The Test 
Center’s activities are concerned with flight and laboratory testing and evaluation 
of guided missiles and their components. 

Each of the establishments has current need for qualified personnel in a variety 
of scientific fields to staff its laboratories. Recently completed at the Naval 
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Ordnance Test Station is Michelson Laboratory at a cost of $6,000,000. Many 
more millions of dollars have been spent in equipment and facilities. Additional 
construction and facilities are planned for both the Air Missile Test Center and 
the Electronics Laboratory. 

The work programs of the laboratories are planned, directed and accomplished 
under the direction of an outstanding staff of civilian scientists. Extensive use 
is made of the council method of operation. Constant liasion is maintained with 
other research organizations, universities, scientific associations, and outstanding 
authorities throughout the nation. 

Professional positions are in the career service of the Federal government under 
Civil Service laws. Examinations are now open in the three scientific establish- 
ments in the following professional fields: Chemist, Mathematician, Metallurgist, 
Meteorologist, Physicist, Statistician, Scientific Research Administrator and 
Scientific Staff Assistant. 

Examinations are also open in the following branches of the Engineering 
profession: Aeronautical, Chemical, Civil, Electrical, Electronics, General, 
Industrial, Material, Mechanical, Metallurgical, Ordnance, Safety and Structural. 

Salaries for most of the positions range from $3397 to $9975 per annum. 
Salaries are predicated on the level of ability, knowledge and experience required 
to effectively discharge the duties of a specific position. 

Further information may be obtained from the Navy Department Joint Board 
of U.S. Civil Service Examiners, 1030 East Green Street, Pasadena 1, California. 


(in eR me 


Reorganization of Philosophy of Science Association 


The Philosophy of Science Association has been reorganized with Philipp Frank 
of Harvard University as President; C. West Churchman of Wayne University, 
Detroit, as Secretary-Treasurer. 

The following are members of the Governing Committee: Gustav Bergmann, 
State University of lowa; Thomas A. Cowan, Wayne University; Clyde Kluck- 
hohn, Harvard University; Sebastian Littauer, Columbia University; F. S. C. 
Northrop, Yale Yniversity. 

The official journal of the Association is the Philosophy of Scicnce of which 
Professor C. West. Churchman is Acting Editor. Manuscripts should be sent to 
the Acting Editor. 

Applications for membership may be sent to the Secretary-Treasurer. Dues 
are $5.00 a year. 

The Association encourages the establishment of local groups in the philosophy 
of science. 

a 
Columbia University Conference on Industrial Experimentation. 


The School of Engineering of Columbia University in the City of New York 
announces an Intersession Five-day Intensive Training Conference on Industrial 
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Experimentation to be offered September 14-18, 1948 by the Department of Ip- 
dustrial Engineering in cooperation with the Department of Mathematical Sta- 
tistics of the Graduate Faculty of Political Science. 

The lecturing will be shared by Professors 8. B. Littauer and J. Wolfwitz and 
a staff of special lecturers drawn from industry. 

A descriptive brochure will be ready for mailing in the latter part of July. 
For further details, interested persons may communicate directly with Professor 


S. B. Littauer, Department of Industrial Engineering, Columbia University, New 
York 27, New York. 


OS ne RI a 


New Members 


The following persons have been elected to membership in the Institute 
(December 1, 1947 to February 28, 1948) 


Angulo, Walter J., B.E. (Johns Hopkins Univ.) Graduate student at Johns Hopkins Uni- 
versity, 5229 Beaufort Ave., Baltimore 15, Maryland. 

Beard, Helen P., Ph.D. (Mass. Institute of Tech.) Assistant Professor of Mathematics, 
Newcomb College, New Orleans 18, Louisiana. 

Blomquist, Nils G., (Univ. of Stockholm) Statistician, Sverige Reinsurance Company, 
Aladdinsvagen 47, Smedslatten, Sweden. 

Bodwell, Charles A., M.S. (Univ. of Michigan) Graduate student at the University of 
Michigan, Box 773, West Lodge, Ypsilanti, Michigan. 

Burnett, Jean, M.S. (Mich. State College) Instructor in Mathematics, Michigan State 
College, 702 Cherry Lane, East Lansing, Michigan. 

Burton, Robert E., Student u: “udichigan University, 1239 Atkinson Avenue, Detroit 2, 
Michigan. 

Byrd, Paul F., M.S. (Univ. of Chicago) Meteorologist, U.S.A.F., Weather Detachment, 
Lockbourne Air Base, Columbus 17, Ohio. 

Cernuschi, Felix, Ph.D. (Univ. of Cambridge) Professor at the University of Montevideo, 
Asociacion Uruguaya de Estadistica, Av. Agraciada 1464, Montevideo, Uruguay. 
Connor, William S., Jr., M.A. (Univ. of North Carolina) Associate Professor of Economics, 

University of Kentucky, College of Commerce, Lexington, Kentucky. 
Dalenius, Tore, Fil. kand. Hastholmsvagen 16, Stockholm, Sweden. 
Davis, Roderic C., M.S. (Calif. Institute of Tech.) P-6 Mathematician, Head of Assessment 
Section, P.O. Box N-467, N.O.T.S., Inyokern, California. 
Dvoretzky, Aryeh, Ph.D. (Hebrew Univ., Jerusalem) Research Fellow at Hebrew Uni- 
versity, Jerusalem, % American Friends Hebrew University, 9 East 89th St., New York. 
Gardner, Robert S., M.S. (Tulane Univ.) Instructor in Statistics, Mathematics Department 
at Ohio State University, 215 W. Eleventh Avenue, Columbus, Ohio. 
Goins, Mary, M.S. (Univ. of Mich.) Assistant Professor of Mathematics, Marshall College, 
435 Ninth Avenue, Huntington, West Virginia. 
Hratz, Joseph A., B.A. (St. Ambrose College) Instructor of Mathematics, St. Ambrose 
College, Davenport, Iowa. 
Kelly, Harriet J., Ph.D. (Univ. of Iowa) Research Associate, Head of Statistical Dept. 
Children’s Fund of Michigan, 660 Frederick, Detroit 2, Michigan. 
Kincaid, Wilfred M., Ph.D. (Brown Univ.) Instructor in Mathematics, University of 
Michigan, Ann Arbor, Michigan. 
Kish, Leslie, B.S. (College of the City of N. Y.) Senior Sampling Statistician, 701 Mt. Pleas- 
ant, Ann Arbor, Michigan. 
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Lehr, Marguerite, Ph.D. (Bryn Mawr, Pa.) Associate Professor of Mathematics, Bryn Mawr 
College, Cartref, Bryn Mawr, Pennsylvania. 

Loizelier, Enrique Blanco, M.A. (Madrid Univ.) Professor of Statistics, Faculty of Fceo- 
nomics, Madrid, University, Nervion 4, Madrid, Spain. 

Lorenzo, Cesar M., M.A. (American Univ., Wash., D. C.) Statistician, Food and Agriculture 
Organization of the United Nations, 1735 De Sales Street, N.W., Washington 6, D. C. 

Lott, Fred W., Jr., M.A. (Univ. of Mich.) Teaching Fellow at the University of Michigan, 

233 Malden Court, Willow Run, Michigan. 

Mantel, Nathan, B.S. (City College of New York) Biostatistician, U.S. Public Health Serv- 
ice, 4421 Lily Ponds Dr., N.E., Washington 19, D.C. 

Marrian, Dixon M., A.M. (Columbia Univ.) Master at the Gilman Country School, Balti- 
more, Md., 1506 Shadyside Road, Baltimore 18, Maryland. 

Martins, Octavio Augusto L., M.A. (Columbia Univ.) Tecnico de Educacao, Department of 
National Education, Rio de Janeiro, Rua Figueiredo de Magalhaes 27, Apt. 702, Rio de 
Janeiro, Brazil. 

Meyer, Herbert Albert, Ph.D. (Univ. of Iowa) Associate Professor of Mathematics, 2016 
West Columbia, Gainesville, Florida. 

Nikitich, Nicholas, B.C.S. (New York Univ.) Timekeeper, 20 Featherbed Lane, New York 
52, New York. 

Oakland, Gail Barker, M.A. (Univ. of Minnesota) Associate Professor of Statistics, Uni- 
versity of Manitoba, Winnipeg, Canada. 

Olkin, Ingram, L.S8. (College of City of N. Y.) Graduate Student at Columbia University, 
245 Fort Washington Ave., New York 32, New York. 












REPORT ON THE NEW YORK MEETING OF THE INSTITUTE 


The thirty-third meeting of the Institute of Mathematical Statistics was held 
at Columbia University, New York City, New York on Wednesday afternoon 
and Thursday, April 14 and 15, 1948. The meeting was attended by 158 persons 
including the following 78 members of the Institute: 











M. Afzal, L. A. Aroian, R. M. Auer, W. D. Baten, R. E. Bechhofer, J. H. Bushey, J. M. 
Cameron, B. H. Camp, G. C. Campbell, 8. D. Canter, Manuel Cynamon, Tore Dalenius, 
J. I. Daly, J. L. Doob, C. W. Dunnett, Aryeh Dvoretzky, Churchill Eisenhart, Benjamin 
Epstein, M. W. Eudey, D. A. Fraser, M. A. Geisler, Mary Goins, H. H. Goode, E, J. 
Gumbel, M. H. Hansen, Mina Haskind, L. H. Herback, 8. M. Ikhtiar-ul-Mulk, Seymour 
Jablon, L. F. Knudsen, Jack Laderman, Howard Levene, 8. B. Littauer, F. M. Lord, 
Irving Lorge, Eugene Lukacs, W. G. Madow, Sophie Marcuse, Robert Mirsky, E. B. 
Mode, D. J. Morrow, Frederick Mosteller, D. N. Nanda, P. M. Neurath, G.*E. Noether, 
M. L. Norden, Ingram Olkin, P. S. Olmstead, A. L. O’Toole, Katharine Pease, E. J. Pit- 
man, W. A. Reynolds, J. S. Rhodes, H. E. Robbins, H. G. Romig, Ernest Rubin, Herman 
Rubin, P. J. Rulon, Frank Saidel, G. R. Seth, M. A. Schlorek, 8S. S. Shrikhande, Rose- 
dith Sitgreaves, Milton Sobel, Emma Spaney, F. F. Stephan, B. R. Suydam, Henry 
’ Teicher, J. W. Tukey, A. Wald, H. M. Walker, J. E. Walsh, 8S. S. Wilks, Dzung-shu Wei, 
Lionel Weiss, Jacob Wolfowitz, C. A. Wright, Mohammad Yusuf. 
















The Wednesday afternoon session, Professor 8. B. Littauer of Columbia Uni- 
versity presiding, was devoted to the following two invited addresses: 















1. Incomplete Block Designs 
Professor R. C. Bose, Caleutta University and the University of North Carolina 

2. Non-Parametric Inference 

Professor J. G. Pitman, University of Tasmania and Columbia University 





The Thursday morning session, Professor Hobart Bushey of Hunter College 
presiding, consisted of a Symposium on Scales of Measurement at which two 
invited papers: 


1. The Development of Psychological Scaling Techniques 
Professor Harold Gulliksen, Princeton University 

2. A Generalized Model for Scales 

Professor Paul Lazarsfeld, Columbia University 


were followed by prepared discussion by Professors Phillip Rulon of Harvard 
University and John Tukey of Princeton University. 

The Thursday afternoon session, Dr. Harry G. Romig of Bell Telephone 
Laboratories presiding, was devoted to the following contributed papers: 















1. Optimum Character of the Sequential Probability Ratio Test 
Professors Abraham Wald and Jacob Wolfowitz, Columbia University 
. Multi-parameter Sequential Estimation 
Mr. G. R. Seth, Columbia University 


to 
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. The Distribution of a Definite Quadratic Form 
Professor Herbert Robbins, University of North Carolina 

. The Moments and Cumulants of the Product of 2, 3, or 4 Dependent Variables (Prelim- 
inary Report) 
Professor Leo A. Aroian, Hunter College 

. Generalization to N Dimensions of Inequalities of the Tchebycheff Type 
Professor Burton H. Camp, Wesleyan University 

. On the Power Function of a Sign Test Formed by Using Subsamples 
Dr. John E. Walsh, Project Rand 

. The Distribution of T?, a Multivariate Generalization of the F-test. 
Miss Dorothy J. Morrow, University of North Carolina. 

. Approximate Confidence Points (Preliminary Report) 
Professor John Tukey, Princeton University 


At all of the sessions there was active discussion from the floor. 
On Wednesday evening members and guests had dinner at the Men’s Faculty 
Club. 


S. B. LirravEr 
Assistant Secretary 


if 
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SANKHYA 


The Indian Journal of Statistics 
Edited by P. C. Mahalanobis 


Vol. VIII, Part 3, 1947 


On some analogues of the amount of information and their use in statistical 

estimation A. BHATTACHARYYA 
The existence of collectives in abstract space... C. F. Kossack 
On recursion formulae, tables and Bessel function populations associated with 

the distribution of Classical D?-Statistic bast oes er P. K. Bose 
On a resolvable series of balanced incomplete block designs... . R. C. Bose 
Notes on testing of composite hypotheses. .. S. N. Roy 
On the general law of demand for raw jute. . P. CHATTERJEE 
Miscellaneous Notes 


Annual subscription : 30 rupees 
Inquiries and orders may be addressed to the 
Editor, Sankhya, Presidency College, Calcutta, India. 
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Journal of the Econometric Society 
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Published Quarterly Subscription to Nonmembers: $9.00 per year 


The Econometric Society is an international society for the advancement of economic theory in ita 
relation to statistics and mathematics. 


Subscriptions to Econometrica and inquiries about the work of the Society and the procedure in 
applying for membership should be addressed to Alfred Cowles, Secretary and Treasurer, The Econ- 
ometric Society, The University of Chicago, Chicago 37, Illinois, U.S.A 
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